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ABSTRACT 

A  novel  approach  to  multiresolution  analysis  based  on  reproducing  kernel  particle  methods 
(RKPM)  and  wavelets  is  presented.  The  concepts  of  reproducing  conditions,  discrete 
convolutions,  and  multiple  scale  analysis  are  described.  By  means  of  a  newly  proposed  semi- 
'  discrete  Fourier  analysis,  RKPM  is  further  elaborated  in  the  frequency  domain,  and  the 
interpolation  estimate  and  the  convergence  of  Galerkin  solutions  are  given.  The  elimination  of  a 
mesh,  combined  with  the  properties  of  the  dilation  and  translation  of  a  window  function, 
multiresolution  analysis,  wavelet-based  error  estimators,  and  edge  detection  brings  about  a  new 
generation  of  hp  adaptive  methods.  In  addition,  this  class  of  multiple  scale  reproducing  kernel 
particle  methods  is  particularly  suitable  for  problems  with  large  deformations,  high  gradients,  and 
high  modal  density.  The  current  application  areas  of  RKPM  include  structural  acoustics,  strucmral 
dynamics,  elastic-plastic  deformation,  computational  fluid  dynamics  and  hyperelasticity. 

1.  INTRODUCTION 

Our  emphasis  in  this  paper  is  on  the  development  of  meshless  methods  for  the  accurate 
prediction  of  the  behavior  of  a  complex  engineering  system  that  involves  a  wide  spectrum  of 
frequencies  and  wave  numbers.  In  particular  this  paper  addresses;  how  does  one  effectively  deal 
with  engineering  systems  characterized  by  multiple  temporal  and  spatial  scales!  In  addition,  can 
accurate  interpolation  functions  be  constructed  with  the  help  of  multiresolution  analysis  concepts 
so  that  the  response  can  be  separated  into  multiple  frequency/wave  number  bands  for  a  better 
representation  of  the  computed  physics!  Our  proposed  approach  is  to  employ  signal  processing 
theories  to  attack  this  difficult  problem. 

Recently,  several  different  meshless  methods  have  been  proposed,  including  Smooth 
Particle  Hydrodynamics  (SPH)  method  [Monaghan  (1988),  Gingold  and  Monaghan  (1977), 
Attaway.  Heinstein,  Mello  and  Swegle  (1993),  Johnson,  Peterson  and  Stryrk  (1993),  Libersky 
and  Petschek(1990)],  Diffuse  Elements  Method  (DEM)  by  Nayroles,  Touzot  and  Villon  (1992), 
Element  Free  Galerkin  (EFG)  [Belytschko,  Lu.  and  Gu  (1994a,b,  1995a),  Belytschko,  Krongauz, 
Fleming.  Organ,  and  Liu  (1995b)  and  Lu,  Belytschko  and  Gu  (1994)],  Particle  In  Cell  methods 
(PIC)  [Sulsky,  Chen,  and  Schreyer  (1992)],  Reproducing  Kernel  Particle  methods  (RKPM)  [Liu, 
Adee  and  Jun  (1993b),  Liu,  Jun,  Li,  Adee  and  Belytschko  (1995a),  Liu,  Jun  and  Zhang  (1995b), 
Liu.  Chen.  Jun.  Chen,  Belytschko,  Pan,  Uras  and  Chang  (1995c),  Liu,  Jun,  Sibling.  Chen  and 
Hao  (1995d).  Liu,  Li.  and  Belytschko  (1995e),  Liu.  Chen  and  Uras  (19950,  Liu  (1995)  and 


Shodja,  Mura  and  Liu  (1995)  and  Liu,  Chen  and  Hao  (1996)],  wavelet  particle  methods  [Liu  and 
Obersti-Brandenburg  (1993a),  Liu  and  Chen  (1995)],  hp  clouds  [Duarte  and  Oden  (1995)], 
partition  of  unity  finite  element  method  [Babuska  and  Melenk  (1995)],  finite  points  method 
[Ohate,  Idelsohn  and  Zienkiewicz  (1995)]  and  free  mesh  method  [Yagawa,  Yamada  and  Kawai 
(1995)].  These  methods  all  have  unique  advantages  and  disadvantages  of  their  own.  On  most  of 
these  methods,  significant  advances  need  to  be  made  before  robust  treatments  of  general  classes  of 
problems  will  be  possible. 

Although  some  of  these  methods  are  currently  being  used  in  the  design  and  analysis  of 
large  scale  engineering  systems,  it  is  stiU  in  infancy  in  answering  the  above  raised  questions.  Even 
though  beautiful  computer  graphics  can  be  generated  to  illustrate  the  numerical  results,  Fourier  type 
analysis  is  still  required  to  study  the  physics  of  the  computed  response  and  to  obtain  the  necessary 
physical  insights.  Moreover,  since  a  Fourier  spectrum  does  not  provide  a  complete  time  or  space 
localization  information,  it  is  desirable  to  have  a  short  time  (or  space)  Fourier  spectrum  available  at 
any  desirable  time  and/or  spatial  location.  However,  all  the  computations  are  performed  in  time 
or/and  spatial  domains.  This  motivates  us  to  employ  a  flexible  time-frequency  and/or  space-wave 
number  window  function  to  construct  global  shape  functions  that  can  separate  the  response  into 
different  scales  (frequencies  and  wave  numbers).  With  this  construction,  Fourier  type  analysis  of 
the  computed  response  is  eliminated,  since  the  frequency/wave  number  bands  (or  scales)  are 
included  in  the  shape  functions. 

Our  intention  here  is  to  introduce  multiple  scale  methods  which  are  based  on  discrete  and 
continuous  reproducing  kernels,  wavelets,  and  integral  window  transforms  and  to  address  some  of 
the  fundamental  issues  related  to  the  development  of  multiple  scale  meshless  methods.  It  is  noted 
that  because  of  the  time/space-frequency/wave  number  localization  via  a  flexible  window  function, 
the  interpolation  functions  are  constructed  so  that  the  response  can  be  separated  into  multiple 
frequency/wave  number  bands  for  a  better  representation  of  the  computed  physics. 

This  class  of  methods  permits  the  response  of  a  system  to  be  separated  into  different  scales, 
with  wave  numbers  corresponding  to  spatial  scales  and/or  frequencies  corresponding  to  temporal 
scales;  so  that  the  response  of  each  scale  can  be  examined  separately.  Through  this  multiresolution 
analysis,  the  physical  interpretation  of  the  computed  results  can  be  immediately  synthesized,  as  in 
classical  Fourier  analysis;  however,  without  the  restrictions  of  the  classical  Fourier  analysis,  which 
is  mostly  restricted  to  simple  geometries  and  discrete  frequencies/wave  numbers. 
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This  development  can  be  likened  to  constructing  a  microscope  with  a  flexible  space-wave 
number  and  a  localized  window  function  which  translates  and  dilates  in  space  and  time  to  cover  the 
entire  domain  of  interest.  This  microscope  can  magnify,  examine,  and  record  the  image  at  vanous 
scales  and  frequencies  locally  within  the  support  of  the  window  function.  The  degree  of 
magnification  will  depend  on  the  power  of  the  microscope,  a  flexible  space-scale  and  time- 
frequency  window  function.  This  characterization  of  the  response  is  performed  through  the  integral 
window  transform.  Localization  can  be  achieved  by  contracting  the  flexible  multiple-scale  window 
function.  The  zoom  in  and  zoom  out  capability  of  the  window  function  is  especially  useful  in 
examining  complex  flow  phenomena,  such  as  flow  induced  vibration,  dynamic  stability  of  flow- 
structure  interaction,  turbulent  structures,  structural  acoustics,  and  high  frequency  structural 
dynamics  response. 

Multiple  scale  methods 

Multiple  scale  analysis  [Liu,  Zhang  and  Ramirez  (1991),  Liu  and  Haeussermann  (1992)] 
has  its  origin  in  signal  analysis.  Wavelet  analysis  [Beklkin,  Coifman,  Daubechies,  Mallat,  Meyer, 
Raphael  and  Ruskai  (1992),  Chiu  (1992),  Mallat  (1989)  and  Strang  (1989),  Williams  and 
Amararnnga  (1994)  and  Daubechies  (1992)]  is  a  contemporary  science  in  image  processing  that  has 
promise  in  computational  mechanics.  However,  one  major  drawback  in  its  application  to 
computational  mechanics  is  its  inability  to  handle  large  deformation  and  complex  boundary 
conditions.  0ne  of  the  key  successes  of  reproducing  kernel  particle  methods  (RKPM)  [Liu,  Jun, 
Li,  Adee  and  Belytschko  (1995a),  Liu,  Jun  and  Zhang  (1995b),  Liu,  Chen,  Jun,  Chen, 
Belytschko,  Pan,  Uras  and  Chang  (1995c),  Liu,  Jun,  Sibling,  Chen  and  Hao  (1995d),  Liu,  Li, 
and  Belytschko  (1995e),  Liu,  Chen  and  Uras  (1995f)  and  Liu,  Chen  and  Hao  (1996)]  is  the 
formulation  of  the  boundary  correction  function  to  scaling  functions  and  wavelets.  Hence,  unlike 
the  usual  wavelet  and  smooth  particle  hydrodynamics  (SPH)  analysis,  no  artificial  boundanes  are 

needed  in  RKPM. 

For  computational  mechanics,  a  discretization  of  the  system  is  inevitable.  When  a  system  is 
discretized,  aliasing  (commonly  refers  to  as  high-frequency  replicas)  is  introduced  into  the 
response.  For  a  complex  system,  aliasing  may  interact  with  the  high  frequency  part  of  the  response 
and  it  becomes  impossible  to  separate  them  clearly  [Liu  and  Chen  (1995)].  In  this  case,  reducing 
the  effect  of  aliasing  is  a  major  step  to  improving  the  solution.  Moreover,  using  this  local  aliasmg 
infonnation.  local  r^inement  or  hp-like  adaptive  refinement  without  a  mesh  can  be  carried  out 
without  the  help  of  the  exact  solution. 
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Liu  and  Chen  (1995)  developed  an  error  estimation  technique  based  on  multiresolution 
analysis,  which  is  especially  useful  for  local  refinement  and  convergence  studies.  The  flexible 
space-scale  window  function  can  be  constructed  to  resemble  the  weU-known  unstructured  multi¬ 
grid  and  hp-adaptive  finite  element  methods.  However,  the  multiple  scale  adaptive  refinements  are 
performed  simply  by  inserting  nodes  into  the  highest  wavelet  scale  solution  region,  and  at  the  same 
time  narrowing  the  window  function.  Hence,  hp-like  adaptive  refinements  can  be  performed 
without  a  mesh.  An  energy  error  ratio  parameter  is  introduced  as  a  measure  of  aliasing  error,  and 
critical  dilation  parameters  are  determined  for  a  class  of  spline  window  functions  to  obtain  optimal 
accuracy.  With  this  parameter,  we  are  able  to  separate  the  numerical  noise  from  the  high 
frequency/wave  number  component  of  the  physical  phenomena.  In  a  traditional  numerical  analysis, 
both  the  numerical  noise  and  the  high  frequency  component  of  the  physics  are  damped  out  by  the 
addition  of  artificial  viscosity. 

Motivation  of  multiresolution  analysis 

A  simple  illustration  of  multiresolution  analysis  is  given  in  Fig.  1.1.  The  left  hand  top 
comer  photo  presents  the  highest  resolution  (256  by  256  pixels).  This  photo  image  is  digitally 
separated  into  a  set  of  consecutive  octave  (power  of  2)  scales  via  a  RKPM  wavelet  transform 
decomposition  (see  Sections  7  and  8).  The  low  and  high  scale  images  of  the  original  photo  are 
presented  below  the  original  photo  on  the  left  hand  side  of  Fig.  1.1.  As  can  be  seen,  the  low  scale 
image  does  .not  contain  any  high  frequency  components,  whereas,  the  high  scale  image  captures 
only  the  sharp  edges  of  the  original  image.  Via  this  multiresolution  analysis,  each  scale  of  the 
image  can  be  studied  separately  and  enhancement  or  modification  can  be  performed  by  simply 
changing  the  local  wavelet  coefficients  at  the  desired  scale  and  location. 


Figure  1.1  Multiresolution  RKPM 
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In  a  classical  Fourier  analysis,  in  which  the  sine  and  cosine  waves  do  not  have  a  compact 
support,  such  modifications  require  the  computation  of  a  new  set  of  Fourier  coefficients  which  can 
be  very  costly.  Presently,  this  is  a  common  practice  in  the  numerical  study  of  the  computed 
physical  response.  In  wavelet  analysis,  unlike  the  classical  Fourier  analysis,  only  a  few  wavelet 
coefficients  need  to  be  computed  because  of  the  space-scale  localization  process.  In  three 
dimensions,  this  localization  process  will  involve  a  space-wave  number  and  time- frequency 
localized  window. 

In  our  research,  we  borrow  this  powerful  concept  of  multiresolution  analysis  to  separate 
the  numerical  noise  from  the  high  scale  components  (frequencies  and  wave  numbers)  of  the 
compressible  flow-structure  phenomena.  Similar  to  the  wavelet  decomposition  of  the  photo  image 
presented  earlier,  the  actual  computed  solution  of  problem  with  a  discontinuity  using  the  multi¬ 
scale  RKPM  hp-like  adaptive  refinement  [Liu  and  Chen  (1995)],  the  low  scale  (scaling  function) 
and  the  high  scale  (wavelet)  parts  of  the  total  solution  are  presented  on  the  right  hand  side  of  Fig. 
1.1.  As  can  be  seen,  the  low  scale  solution  consists  of  the  low  wave  number  solution,  whereas, 
the  high  scale  solution  consists  of  a  mixture  of  the  high  wave  number  approximation  of  the  physics 
and  the  undesirable  numerical  noise,  or  the  aliasing.  Analogous  to  the  low  scale  photo  image,  the 
low  scale  solution  depicts  only  the  smooth  part  of  the  solution  and  the  edges  of  the  discontinuity 
(represented  by  a  high  wave  number  wavelet  solution)  have  been  removed.  However,  the  edges  of 
the  discontinuity  can  be  located  via  this  high  scale  component.  This  shows  that  not  only  the 
location  of  the  discontinuity  can  be  detected,  numerical  noise  can  also  be  minimized  by  aliasing 

control. 

Section  2  describes  the  scaling  functions  and  wavelets  associated  with  multiresolution 
analysis.  Two  different  methods  for  constructing  the  Daubechies  scaling  function  are  reviewed 
[Daubechies  (1992)].  In  Section  3,  the  orthogonality  conditions  for  the  scaling  functions  and 
wavelets  are  studied  in  the  Fourier  domain.  A  systematic  procedure  of  constructing  orthogonal 
scaling  functions  and  orthogonal  wavelet  functions  based  on  the  Fourier  transform  is  also  derived. 
The  presentations  in  Sections  2  and  3  hopefuUy  provide  an  easier  exposition  of  the  construction  of 
the  scaling  functions  and  wavelets. 

The  concepts  of  scaling  functions  and  wavelets  are  then  used  to  construct  reproducing 
kernel  methods  (RKM).  The  discrete  form  of  a  reproducing  kernel  method  is  called  a  reproducing 
kernel  particle  method  (RKPM)  and  forms  the  basis  of  applying  these  concepts  in  computational 
mechanics.  In  Section  4,  the  scaling  function  of  RKM  is  constructed  by  relaxing  the  orthogonality 
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condition.  Consequently,  a  nonuniform  discretization  of  the  reproducing  kernel  can  be  achieved  via 
a  particle  method.  The  concepts  involved  in  the  development  of  RKM  from  Taylor  series,  using 
multi-index  notation  are  introduced  in  Section  4.  An  orthogonal  window  function  for  RKM  is  also 
briefly  proposed  in  the  same  section.  The  formulation  of  RKPM,  the  discrete  form  of  RKM,  is 
presented  in  Section  5. 

The  study  of  RKM  and  RKPM  formulation  in  the  Fourier  space  is  discussed  in  Section  6. 
Interpolant  estimates  and  convergence,  and  their  relationship  with  the  wavelet  solution  are  explored 
in  Section  7.  Multiresolution  analysis  for  RKPM  is  reviewed,  and  the  application  of  RKPM  in 
edge  detection  and  adaptivity  in  computational  mechanics  are  outlined  in  Section  8.  The 
applications  of  RKPM  in  large  deformation,  computational  fluid  dynamics  and  strucmral  acoustics 
are  presented  in  Section  9,  followed  by  conclusions. 


2.  SCALING  FUNCTION  AND  WAVELETS 

2.1  Preliminaries 

Multiresolution  Analysis  and  Scaling  Function 

A  multiresolution  analysis  makes  use  of  a  nested  sequence  of  closed  subspace  { Vj}j^2 

/ 

a  function  <p  with  the  following  properties 

(a).  ...  c  V.i  c  VqC  Vi  c--' cl2(R) 


(b) .  u  Vj  =  lHR) 

ieZ 

(c) .  n  V,=  {01 

jeZ 

(d) .  Each  subspace  is  related  to  its  next  finer  level  subspace  by 

0(.x)  €  Vj<^(p{2x)^  Vy+i,  j^Z 

(e) .  The  ttanslates  of  <p{x)  span  the  same  subspace  so 

(^x)  e  Vo  <p{x  -  n)  e  Vq,  for  all  ne  Z 
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For  an  orthogonal  wavelet  basis,  we  have  the  additional  condition: 

(f).  The  function  0(x  -  n)  is  a  mother  function  for  an  orthonormal  basis  of  Vq,  where 


0j„  =  2’>^(tK2jx-n), 

for  all  i,n€Z,  Z  =  {•  •  •,  -2,  -1,  0,  1,  2,  ■  •  •}.  Conditions  (d)  and  (f)  imply  {(pj.n  \neZ}  is  an 
orthonormal  basis  of  Vj  for  all  /  e  Z  . 

Since  0  e  Vq  c  Vi  ^  and  {0i,„ }  forms  an  orthonormal  basis  in  Vi ,  any  function  in  Vo  can 
be  expressed  in  terms  of  the  basis  functions  of  Vi .  Thus, 


oo 

(tKx)=  X  a„(l)(2x-n) 

n  =  -» 


(2.1.2) 


where  a„  are  given  values  for  each  particular  function  <^x),  which  is  often  referred  to  as  the 
“scaling  function”  of  the  multiresolution  analysis. 


Wavelet 

A  new  subspace  Wj.\  associated  with  the  multiresolution  analysis  is  defined  as. 

/ 

Vj=Vj.x®Wj.x  ' 

where  ©  is  a  direct  sum,  such  that  Wj.x  the  orthonormal  complement  of  V,.i  in  Vj,  for  every 
/  e Z.  It  follows  that  for,;  >  /, 


(2.1.4) 


where  all  subspaces  Wj  are  orthogonal  and  that 

©  Wj  =  lHR) 

lez 


(2.1.5) 
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Basically,  for  multiresolution  analysis,  if  a  sequence  of  subspaces  Vj  satisfies  the  properties  (a)  - 
(0,  there  exists  an  orthonormal  wavelet  basis, 

\f/j  „  =  2j'^\lK2jx-  n), 

such  that ,  for  fixed  j,  ( Vj.n  \n€Z}  constitutes  an  orthonormal  basis  for  Wj.  By  virtue  of  (b),  (c) 
and  Eq.  (2.1.5),  this  also  implies  that  the  whole  collection  ( !/(/,„  «  e  Z)  is  an  orthonormal  basis 

for  Furthermore,  from  Eq.  (2.1.3),  the  wavelet  function  yKx)  can  be  expressed  in  terms  ot 

the  scaling  function  0(x)  at  the  next  finer  scale. 


\lKx)  =  bn  (p{2x  -  n) 


Illustration  of  Multiresolution  Analysis  by  the  Hoar  Function 

As  an  example  of  a  scaling  function  and  associated  wavelet,  consider  the  Haar  basis.  The 

Haar  function,  also  known  as  the  box  function,  is  given  as: 

j  1,  0<x<  1  (2.1.7) 

I  0,  otherwise 

I 

From  Figure  2.1a,  it  is  noted  that  the  box  function  satisfies  Eq.  (2.1.2)  with ‘coefficients 
aQ  =  a\  =  \  thus 

(p{x)  =  <p{2x)  +  <p{2x  -  1)  (2.1.8) 

Furthermore,  as  shown  in  Figure  2.1b,  the  Haar  wavelet  satisfies  Eq.  (2.1.6)  with  coefficients 
bQ  =  \  and  foi  =  -1,  so  that: 


VA.t)  =  0(2x)  -  0{2x  -  1) 


(2.1.9) 
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Figure  2.1  Haar  Scaling  Function  and  Wavelet 


2.2  Daubechies  Scaling  Function 

In  order  to  determine  <p{x)  at  a  set  of  points,  the  filter  coefficients  in  Eq.  (2.1.2)  which 
satisfy  conditions  (a)  through  (e)  are  required.  These  conditions  are  derived  from  the  properties  of 
the  scaling  function  which  forms  an  orthogonal  basis. 

Conditions  for  the  Determination  of  Filter  Coejficients 

(i)  The  area  under  the  scaling  function  is  normalized  to  unity  to  ensure  that  each  scaling 

function  of  a  given  shape  is  uniquely  defined.  Thus, 
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Thus,  the  first  condition  on  these  A^-filter  coefficients  is  given  as: 


N-\ 

an  =  2 


n=0 


(2.2.3) 


(ii)  The  scaling  function  has  to  be  orthogonal  to  its  integer  translates,  that  is. 


^  0{x)  0ix  +  k)  dx  =  5ojt,  k  —  0,  1 , ...,  A72- 1  (2. -.4) 

where  Sok  is  the  Kronecker  Delta;  in  the  orthogonality  condition,  N  is  restricted  to  be  an  even 
integer.  Substituting  Eq.  (2.1.2)  into  Eq.  (2.2.4); 

{  ^  a„<l>{2x  -  n)  ^  an<t>{2x  +  2k -n)  dx  =  5ok,  ^  =  0,  1, ...,  A//2-1  (2.2.5) 

n  =  0  «=  0 


and  rearranging  the  indices  for  n=m  +  2k  yield: 


V-l  N-\-2k 

^  fl„</>(2x  -  n)  ^  am+2k<t>i2x  -  m)  dx  =  5oh  k  =  0,  1, 

n  =  0  m  =  -2/c 


M2-1  (2.2.6) 


From  Eq.  (2.2.4),  Eq.  (2.2.6)  results  in; 


NA  N-\-2k 
^  an  ^  am+2k  ^nm 

n  =  0  m  =  -2k 


=  25ok,  k  =  0,\ . M2-1 


(2.2.7) 


Therefore,  the  orthogonality  condition  in  Eq.  (2.2.7)  provides  N/2  equations: 
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(2.2.8) 


Up  to  this  point,  we  have  NH  +  1  equations  to  define  N  coefficients,  so  another  W2  -  1 
equations  are  necessary  to  determine  a  unique  set  of  filter  coefficients.  The  approach  used  is  to 
require  the  scaling  function  to  be  able  to  represent  polynomials  of  order  up  ioN/2-\  exactly.  This 
constraint  leads  to  the  compactly  supported  wavelets  [Daubechies  (1992)]. 

(iii)  Consider  a  given  polynomial  function  with  order  N/2\ 


f{x)  =00+  (X\x  +  02^2  +  ...  +  OAr/2-lX 


NUA 


(2.2.9) 


This  function  can  also  be  represented  by  the  following  expansion: 


N-\ 

f(x)=  X  c„<p{x-n) 

n  =  0 


(2.2.10) 


The  wavelet  can  be  used  to  obtain  conditions  on  the  filter  coefficients  by  taking  the  inner  product 
of  Eq.  (2.2/10)  with  yKx),  thus: 


N-l  , 

(fix),  vKx))  =  X  W ® 

n  =  0 


(2.2.11) 


This  is  because  0{x  -  n)  is  required  to  be  orthogonal  to  the  wavelet  function,  XiKx).  Hence,  using 
Eq.  (2.2.9)  in  the  left  hand  side  of  Eq.  (2.2.1 1)  gives: 


CCoj 


[ 


\lA,x)dx+ai\  x\i/{x)dx  +  ...  +  aM/2-\\  x^^^-^\iKx)dx -0  (2.2.12) 


Eq.(2.2.12)  holds  for  arbitrary  Oa  so  that  the  additional  N/2  conditions  are: 


i; 


x^y/(.x)dx  =  0,  =  0,  1,  2 . N/2  -  1 


(2.2.13) 
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To  relate  Eq.  (2.2.13)  with  the  filter  coefficients  an,  the  Fourier  domain  is  used.  The  Fourier 
transform  of  (pix),  gives: 


f[<l>ix)]  =  0{O})  = 


a„e 


(2.2.14) 


Eq.  (2.2.14)  can  be  expressed  as: 


y=i 


-  CO  TT  CO  ^ 
=  11^  2^ 


;=i 


(2.2.15) 


where-  Pico)  =  For  an  N-filter,  a  wavelet  function  can  be  expressed  as  [Strang 


(1989)]: 


1 

¥ix)=  ^(-l)”^i-n<A(2^^”) 


(2.2.16) 


n=2-N 


Then  the  Foprier  transform  of  Eq.  (2.2.16)  gives: 
¥i(0)  =  ^  ^i-l)''a,.„e‘''‘^^^i^) 

n=2-N 

Recall  the  zero  moment  in  Eq.  (2.2.13), 


c 


\j/ix)cix  =  ^(0)  =  0 


(2.2.17) 


(2.2.18) 


When  0)  =  0.  by  using  Eq.  (2.2.18),  Eq.  (2.2.17)  yields: 


n=2-N 


n=2-N 


(2.2.19) 
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Note  that  0(0)  =  1,  which  can  be  easily  verified  from  the  Fourier  transform  of  Eq.  (2.2.1).  From 
the  definition  of  P{(0),  we  observe  that  P{(0  =  n)  =  ^^(-l)”a„ .  Subsututing  this  relation  into 

Eq.  (2.2.19)  leads  to  P(;r)  =  0  since  ^(0)  =  0.  For  the  ^  moment  in  Eq.  (2.2.13),  it  yields; 


J  x^\lf{x)dx  =  {if  vr^^^(O)  -  0 


(2.2.20) 


The  kth  derivative  of  Eq.  (2.2. 17)  at  ©  -  0  gives: 


(2.2.21) 


n=2-N 


From 


I 

Eq.  (2.2.21)  and  the  kih  derivative  of  P{co),  ,  it  shows  that. 


P^**(;r)  =  0 

Thus,  the  third  condition  in  terms  of  the  a„  can  be  obtained  from  Eq.  (2.2.22). 


(2.2.22) 


(2.2.23) 


It  is  noted  that  the  equation  with  it  =  0  in  Eq.  (2.2.23)  is  redundant  because  it  can  also  be 
found  from  the  combination  of  Eqs.  (2.2.3)  and  (2.2.8).  Thus,  by  excluding  this  redundant 
equation,  Eq.  (2.2.23)  provides  only  N/2  -  1  new  equations.  The  scaling  function  and  wavelet  are 
uniquely  defmed  by  the  coefficients  a„  which  are  obtained  by  solving  the  system  of  (V  equations. 


An  example  of  filter  coefficient  for  N  —  4,  the  Daubechies  D4  scaling  function 


Oq  +  (71*+  <32  ^3  ~  “ 

■>  7  ■>  2 

Aq  +  fl]  +  ^2  ^3  — 


(from  Eq.  (2.2.3)) 

(from  Eq.  (2.2.8)  k  =  0) 


(2.2.24a) 

(2.2.24b) 
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(from  Eq.  (2.2.8)  k=l) 

(2.2.24c) 

^0  ”  +  ^2  ”  ^3  "  ^ 

(from  Eq.  (2.2.23)  ^  =  0) 

(2.2.24d) 

+  2^2  ~  3a3  =  0 

(from  Eq.  (2.2.23)  k=l) 

(2.2.24e) 

In  this  example,  the  redundant  equation  Eq.  (2.2.24d)  can  be  obtained  by. 

-^2  *  Eq.{2.1.2Ab)  +  4  *  Eq.{2.2.24c)  -  Eq.i2.2. 24a)^  (2.2.25) 

The  filter  coefficients,  can  easily  be  obtained  by  replacing  the  nonlinear  equation  00^2  +  =  0 

with  the  linear  equation  Oq -fli +^2 -^3  =0’  which  is  the  redundant  equation.  The  filter 

coefficients  obtained  from  Eq.  (2.2.24a-e)  are  also  known  as  the  coefficients  for  the  Daubechies 
D4  scaling  function.  Table  2.1  lists  the  filter  coefficients  for  the  Daubechies  D4 


n 

an 

0 

0.68301270189222 

1 

1.18301270189222 

2 

0.31698729810778 

3 

-0.188301270189222 

Table  2.1  Filter  coefficients  for  the  Daubechies  scaling  function  D4 

2.3  Construction  of  Scaling  Function 

Given  a  set  of  filter  coefficients  as  defined  in  Section  2.2,  the  corresponding  scaling 
function  can  be  constructed  by  the  iteration  method  or  by  the  recursion  method. 

2.3.1  The  Iteration  method 

The  iteration  method  starts  from  Eq.  (2.1.2); 

00 

0i{x)=  X  On  -  n)  (2.3.1) 


«=  -00 


ri,  0  <  x  <  1 

with  the  box  function  as  (poix),  where  <j>o{x)  -  otherwise 
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Example  1 :  aQ  =  2 

For  ao  =  2,  the  iteration  equation  obtained  from  Eq.  (2.1.2)  yields: 

(^,(x)  =  2  0m(2x)  for  />! 

The  Delta  function  constructed  by  iteration  method  with  aQ=2  is  depicted  in  Fig2.2.  It  is 
obtained  when  /  ^ 


Figure  2.2  The  Delta  function  constructed  by  iteration  with  ao-2 


Example  2  t  ao  =  ^,  a\  -  I,  a2  -  2 

Por  =  1  ^2^  =  1 ,  a2  =  the  iteration  equation  is  given  as: 

m  =  i  (l>i.i(2x  )  +  </)/.i(2x-l  )  +  ^  <^,-i(2x-2  )  for  /  >  1 

The  hat  function  constructed  by  iteration  method  with  00  =  ^,01  =  1,02  =  ^  is  shown  in  Fig.  2.3 
It  is  obtained  when  1  — >  . 
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Example  3:  D4  filter  coefficients 

The  coefficients  in  Eq.  (2.2.2 1 ),  lead  to  the  D4  scaling  function.  The  plot  of  the  D4  scaling 
function  as  constructed  by  iteration  is  shown  in  Fig.  2.4b. 

Remark;  The  iteration  method  is  simple  and  straight  forward,  although,  it  is  not  computationally 
efficient. 


2.3.2  The  Recursion  Method 

The  crux  of  the  recursion  method  hinges  on  the  fact  that  the  scaling  function  (/)(x)  is  known 
at  integer  points,  x  =  j.  Thus,  recursion  gives  the  values  of  (p(x)  at  the  half  integer  pomts  and  at  the 

quarter  integer  points  for  the  next  recursion.  Repeating  this  process  will  yield  all  dyadic  points 
{iKj  ^  0).  This  leads  to  a  fast  algorithm  for  wavelet  calculations.  To  compute  (p{x)  at  the  integer 

points,  consider  the  expansion  of  Eq.  (2.1.2): 

0(x)  =  aQ(l)(2x)  +  ai<t>{2x-  1)  +  ...  +  .  \0{2x -N  ■¥  1)  (2.3.2) 


For  N  point  values  of  0{x),  Eq.  (2.3.2)  can  be  written  in  matrix  form: 


^0 

0 

0  ••• 

0 

0 

0 

■  0iO) 

0(0) 

/ 

Oq  0 

0 

0 

0{l) 

0(1) 

04 

03 

02 

Oq 

0 

0(2) 

0(2)  • 

^6 

as 

04  ••• 

0 

0i^) 

— 

0(3) 

0 

0 

0  ••• 

%-2 

aN-3 

0{N-2) 

0(N-2) 

0 

0 

0  ••• 

0 

0 

i 

1 

JiN-l)_ 

(M-I)<I>  =  0 


The  vector  values  of  O  at  the  integer  points  is  the  eigenvector  of  M  corresponding  to  the 
eigenvalue  1.  As  in  all  eigenvalue  problems,  this  vector  is  not  uniquely  defined.  An  additional 
normaUzing  condition  obtained  from  Eq.  (2.2.1)  is  used  to  determine  a  unique  eigenvector.  Thus. 


^  0{n)  =  l,  neZ 


(2.3.5) 
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Once  we  know  the  values  of  <p{x)  at  the  integer  points,  then  the  values  of  (p(x)  at  the  half  integer 
points  can  be  computed  by  Eq.  (2.1.2); 

n  -  -oo 

Example  4:  D4  scaling  function 

As  an  example  the  recursion  method,  the  values  of  D4  at  four  integer  points  are  calculated 
by  Eq  (2.3.3): 


0 

0 

o' 

■0(0)' 

■0(0)' 

02 

% 

0 

0(1) 

0(1) 

0 

02 

^1 

0(2) 

0(2) 

_0 

0 

0 

^3. 

.0(3). 

The  results  are: 

,^0)  =  0,  0(1)=J-^,  0(2)  0(3)  =0 

Then,  from  Eq.  (2.3.6),  the  values  of  <p{x)  at  the  half  integer  points  are  given: 

0(^)=2^.  0(|)  =  O.  0(|)  =  ^^ 

A  comparison  of  the  D4  scaling  functions  obtained  by  iteration  and  recursion  are  shown  in  Figure 
2.4.  Apparently,  the  construction  from  recursion  is  more  efficient  and  accurate  than  from  iteration. 


(a) 
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Figure  2.4  The  D4  scaling  function  constructed  by  recursion  (a)  and  by  iteration  (b) 


Derivative  of  Scaling  Functions 

The  derivative  of  the  scaling  function  can  be  computed  by  differentiating  Eq.(2.1.2)  with 
respect  to  jc,  thus. 


oo 

<t>'{x)=  '^2a„(l>'{2x-n) 

nzz—oo 


(2.3.7) 


This  leads  to  the  same  eigenvalue  problem  as  for  the  scaling  function,  only  now  we  are  looking  for 
the  eigenvector  corresponding  to  the  eigenvalue  1/2  (c.f.  Eq.  (2.3.3)) 


'flo 

0 

0 

O' 

>'(0)' 

>'(0)' 

aj 

4o 

0 

</>'(!) 

</>'(!) 

0 

43 

(p'a) 

0'(2) 

_0 

0 

0 

fl3_ 

yo)_ 

yo)_ 

The  normalizing  condition  for  the  derivative  arises  from  the  reproducing  property  of  the  scaling 
function.  That  is; 


y  (pix  -  y)  dy  =  x 


(2.3.9) 


Differentiating  Eq.  (2.3.9)  with  respect  to  a:,  yields: 


20 


j  y  (p  (x  -  y)  dy  =  1 

Thus, 

(l>'{x-xj)  =  l  (2.3.11) 

©o 

Eq.  (2.3.1 1)  gives  the  normalizing  condition  for  the  derivative.  The  derivative  of  D4  is  shown  in 
Figure  2.5. 


Figure  2.5  The  derivative  of  D4  scaling  function 


Remark: 

Scaling  functions  which  are  constructed  by  iterative  or  recursion  methods  are  defined 
only  at  dyadic  points  (2A;  >  0).  Therefore,  these  scaling  functions  and  their  derivatives  are  only 
weakly  continuous  and  they  are  often  known  as  fractal  interpolation  functions.  When  y  they 

converge  to  a  continuous  function. 


3.  FOURIER  ANALYSIS  OF  ORTHOGONAL  SCALING  FUNCTIONS 
3.1  Orthogonal  Scaling  Function 

•An  orthoaonai  scaling  function  can  also  be  constructed  by  the  Fourier  transform.  Rewriting 
Eq.  (2.1.2)  as: 
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(}>{x)  =  ^  Crt  (Pl.ni^)  -  ^  0(2x  -  tt)  n  €  Z 


(3.1.1) 


where  c„  are  given  for  each  particular  scaling  function.  The  Fourier  transform  of  0ix)  is  defined  as. 


■I 


(^(ty)  =  <Kx)  dx 


(3.1.2) 


By  changing  variables,  y  =  2x-n  and  x  =  (y+n)/2,  the  Fourier  transform  of  (K2x  -  n)  denoted  by 
F[  ]  is: 

Fm2x  -  n)]  =  I  ^)  dy  =  j  ^Xy)  0(^)  (3.1.3) 

Equation  (3.1.3)  depicts  the  propeny  of  the  Fourier  transform  of  a  "shifted  function  <}>(2x  -  n). 
Thus,  the  Fourier  transform  of  Eq.  (3.1.1)  is: ' 


0(<u)  =  :^'Zc„  (tXf)  e-^  =  md^)  (tXf) 


(3.1.4) 


nr 


where 


V2  n 


(3.1.5) 


Note  that  the  function  mdcS)  represents  the  connection  between  two  different  scales  of  the  scaling 
function  in  the  Fourier  space. 

The  Parseval  identity  and  the  orthogonality  condition  for  the  scaling  function  <f(x)  give  the 
Kronecker- delta  property: 


j*  </)(x)  Cfx-k)  dx=  :^j  <p{0))  F[<pix-k)]  dco  =  5 ok  ’.  keZ 


(3.1.6) 


where  0  denotes  the  complex  conjugate  of  0.  Using  the  property  of  the  Founer  transform  of  a 
shifted  function,  the  left  hand  side  of  Eq.  (3.1.6)  becomes. 


0(x)  0{x-k)  dx  =  :^  0(co)  0(0))  dec 


0(o})0(o})  e‘^d(jt)  = |<^(o>f2;r/)|  je'^do) 


(3.1.7) 


In  deriving  Eq.  (3.1.7),  the  integral  is  subdivided  into  intervals  [2;r/,  2;r(/+l)]  for  /e  Z.  In 
obtaining  Eq.  (3.1.7),  the  periodicity  condition  =  i  has  been  employed.  To  satisfy  the 
condition  in  Eq.  (3.1.6),  we  need 


^  0(ru+2;r/)l  =  1 


(3.1.8) 


which 


1 

is  the  orthogonality  condition  for  a  general  scaling  function  0(x).  Since  \0(o>¥27vl 


')| 


for  arbitrary  0,  the  orthogonal  scaling  function  0  (co)  is  defined  by; 


0  (0))  = 


(3.1.9) 


^  </)(u>f2;r/)l 


Note  that  the  normalization  factor  in  Eq.  (3.1.9)  can  be  expressed  in  terms  of  a  Fourier  series  such 


\0{O)  +  2kI)\ 


_  -inco  ^ 

a„  e  ,  a 


•-tJ  ^ 

271  Jo  , 


r 

\\ko)  +  27:l)f  e‘"^dO)  (3.1.10a,b) 


Thus  the  orthogonal  scaling  function  Eq.  (3.1.9)  becomes: 
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0  (fi))  =  X  ccne-'^^^fpico) 

n 


(3.1.11) 


and  since  the  right  hand  side  of  the  above  is  a  Fourier  transform  of  a  shifted  function,  Eq.  (3. 1.11) 
yields: 


0*(^)  =  X 

n 


(3.1.12) 


3.2  Orthogonal  Wavelet 

* 

In  order  for  v'*(jc)  to  be  an  orthogonal  wavelet  for  an  orthogonal  scaling  function  </»  (.r),  the 
inner  product  of  and  (l)*ix-k)  must  vanish.  That  is 


{\l/\(t>k)  =  <t>  ix-k)  dx  =  0 


(3.2.1) 


Similar  to  Eq.  (3.1.1),  ¥*ix)  can  be  expressed  as: 


¥ 


*ix)  =  X  fn  ihnix)  =  X  ^  -  n) 


(3.2.2) 


where  the  fn  coefficients  are  to  be  determined.  The  Fourier  transform  of  Eq.  (3.2.2)  is. 


¥ico)=^J,f:  <t>  =  (p  {f 


(3.2.3) 


V2 


where 

»>/«») I 

V2  fi 

Using  Parseval's  identify  and  the  property  of  the  Fourier  transform  of  a  shifted  function  Eq. 
(3.1.8).  the  Fourier  transform  of  Eq.  (3.2.1)  is  given  as: 


i//  (co)  (p  (CO)  e‘^‘^dco  =  0 


(3.2.5) 


The  domain  of  the  above  integral  is  subdivided  into  intervals  [Inl,  27t{l+\)]  and  the  integrals  over 
these  intervals  summed.  The  left  hand  side  of  Eq.  (3.2.5)  becomes. 


V  -1-  y/  {co)  0  (co) 
^  27t 


(3.2.6) 


'  2ji 


By  substituting  a  =  a -27d,  Eq.  (3.2.6)  yields: 


J-Y  I  \f/*{co+2Kl)  0  iQ)+2Kl) 
2k ' 


(3.2.7) 


Finally,  using  the  periodicity  of  i.e.  =  l  and  setting  (o-  co,  orthogonality  gives. 


eicok'^  y/  (q)+2kI)  0  {GH-2Kr)do)  =  0 


Jo 

Hence 


(3.2.8) 


^  y/*iccH-2Kl)  0  io>^2K[)  =  0 


(3.2.9) 


Substituting  Eq.  (3.2.3)  and  Eq.  (3.1.4)  into  Eq.  (3.2.9)  yields: 


X  nicif+Kl) 


0  (^^/)1  =0 


(3.2.10) 


From  Eqs.  (3.1.5)  and  (3.2.4),  both  mfico)  and  mdco)  are  periodical  over  2n.  That  is: 
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Regrouping  the  sums  for  odd  and  even  /,  and  using  the  orthogonality  condition  Eq.  (3.1.10),  and 
Eq.  (3.2.11),  Eq.  (3.2.10)  leads  to: 

mf{0)')mc{Q)')  +mf{Q)'+r:)mc{co+K)=0  .  where  0)=^  (3.2.12) 


To  solve  Eq.  (3.2.12),  let  co=  O)  and  assume  m^O))  =  A(<i))  mc{Ct>¥n).  Eq.  (3.2.12)  becomes: 

A(G))  mJBj  mc{0}¥7C)  +  Mo}+n)mc{Qy^rt)  m^ioi+ln)  =  0  (3.2.13) 

From  the  definition  of  mdo)),  Eq.  (3.1.5),  we  have  m,(G))  =  mdoy^ln)  and  mdo))  ^  0, 

so  that  Eq.  (3.2.13)  can  be  simplified  as: 


A(g))  +  l{c»¥n)  =  0 
The  solution  to  Eq.  (3.2.14)  is: 

A(m)  = 
therefore 

/ 

( 

mf{(0)  =  e‘^  mdoy^rr) 


(3.2.14) 


(3.2.15) 


(3.2.16) 


Employing  Eq.  (3.2.16)  in  Eq.  (3.2.3),  the  Fourier  transform  of  an  orthogonal  wavelet  is: 


✓V  * 

y/* {CO)  =  mj{^)  (p  (^)  = 


md^n)  (p  (^) 


(3.2.17) 


Expanding  md^rt)  via  Eq.  (3.1.5)  yields: 


n  n 


^ ±. y.,c; (.1)" 7(f) cli  (-I)'"''  0 

VI  n  2  V2  „ 


(3.2.18) 
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With  Eq.  (3.1.3),  the  inverse  Fourier  transform  of  Eq.  (3.2.18)  gives: 


=  \TX  Cl  <p\2x-n) 

n 


(3.2.19) 


Once  the  coefficients  c*n  are  determined,  the  orthogonal  wavelet  /(x)  can  also  be  found  by  Eq. 
(3.2.19).  Equating  Eq.  (3.2.2)  to  (3.2.19)  gives: 


3.3  Derivation  of  Orthogonal  Wavelets  from  a  Scaling  Function 

In  order  to  derive  an  orthogonal  wavelet  from  a  scaling  function,  Eq.  (3.1.9)  is  rewritten  as: 


(t>  (ty)  = 


(Koi) 


vr  |0(fi)+2;r/)l 


mcif) 


«(f+2«0 


«f) 


VllJ(w2«of  Yl|«(f+2;r;) 


=  (f) 

2  II 


(3.3.1) 


In  deriving  Eq.  (3.3.1).  the  expansion  Eq.  (3.1.4)  is  used  to  represent  «(tB)  and  the  last  term  in 
Eq.  (3.3.1)  can  be  expressed  as  J’lf)  by  “sinS  E<)-  1  S)-  Hence  a  new  ml  for  the  orthogonal 

scaling  function  (p  (co)  can  be  defined  as: 


=  nicim 


<l>{^2nl) 


(3.3.2) 


Substituting  Eq.  (3.1.5)  into  Eq.  (3.3.2)  yield  the  computational  formula  for 


27 


fl 

1 

0(f+2;r/)|' 

2 

0{oyr2Kl) 

i 


(3.3.3) 


It  is  noted  that  for  each  scaling  function  (p{x),  the  coefficients  c„  are  given.  To  represent  m,*  (|)  via 
an  expression  similar  to  Eq.  (3.1.5),  a  new  set  of  coefficient  dn  is  defined  through: 


Z  yZ  n 

Once  me  (^  is  given  by  Eq.  (3.3.3),  the  Fourier  coefficients  d„  appearing  in  Eq.  (3.3.4)  can  be 
computed  as: 


d  =  — 
In 


i: 


fl  mciCD)  dco  = 


VI ^ 


(3.3.5) 


With  d„  determined  and  using  the  expansion  Eq.  (3.1.1),  for  a  given  scaling  function  <p{x),  its 
orthogonal  Scaling  function  is: 


0\x)  =  VJI  d„  0  {Ix-n) 

n 


(3.3.6) 


the  corresponding  orthogonal  wavelet  is  (c.f  Eq.  3.2.19): 

V/*(,r)  =  VlX  ^-”-1  ^*(2x-n)  (3.3.7) 

n 

Substituting  Eq.  (3.1.16)  into  Eq.  (3.3.7),  the  relation  between  an  orthogonal  wavelet  and  a 
scaling  function  that  satisfies  Eq.  (3.1.1)  can  be  expressed  as: 

V/*(.v)  =  ViX  (-!)"■'  ^-"-1  X  00-x-m-n)  (3.3.8) 

n 


3.4  Construction  Procedure 
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The  following  procedures  are  used  to  derive  an  orthogonal  scaling  function  ^  and  an 
orthogonal  wavelet  Xj/*  from  a  scaling  function  0. 


STEP  1: 


Choose  0  so  that: 

1.  0(x)  and  0(0))  to  have  a  reasonable  decay  and  a  finite  support. 

2.  Eq.  (3.1.1)  is  satisfied,  that  is  c„  are  given. 

3.  f  0(x)  dx^O 


(3.4.1a) 

(3.4.1b) 

(3.4.1c) 


STEP  2:  Compute  X  k  ^2  ;r/)| 

i 

The  coefficient,  X  \0(o>^2K[f ,  can  be  expressed  in  terms  of  a  Fourier  series  such  that: 


X  \0(oy^27d'\  =  X  ; 


_=  tL  X  \0(o>^27:i\  do) 

JLTC  I 


(3.4.2a,b) 


Note  that  the  expression  in  Eq.  (3.4.2a,b)  is  essentially  the  same  as  the  last  expression  in  Eq. 
(3.1.7).  Therefore,  Eq.  (3.4.2a,b)  is  simplified: 


hn  -  I  0{x)  0(x-k)  dx 


(3.4.3) 


In  general,  at  least  (m+l)/2  Gauss  points  are  needed  for  the  numerical  integration  of  Eq.  (3.4.3), 
where  m  is  the  highest  order  of  the  product  0(.x)  0{x-k)  in  Eq.  (3.4.3). 

STEP  3:  Compute  — ========  and  coefficients  a„ 

/>/  X  0(6>^-2^/) 


Recall  from  Eq.  (3.1.10a,b) 
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It  is  noted  that 


Vi 


:</)(a>f2;r/)l 


has  been  expressed  in  terms  of  bn  via  Eq.  (3.4.2a,b). 


STEP  4:  Construction  of  the  orthogonal  scaling  function  (p  (x) 
From  Eq.  (3.1.12); 


n 

The  an  coefficients  are  in  general  nonzero  in  Eq.  (3.1.10b).  In  order  to  keep  only  a  few  terms  an, 
we  choose  rapidly  decaying  scaling  functions  so  that  only  a  finite  number  of  terms  is  necessary.  In 
other  words  a„  decays  rapidly  away  from  /i  =  0.  Thus,  we  can  neglect  the  high  order  coefficients 


STEP  5:  Compute  mc(y) 

We  employ  the  scaling  function  (p{x)  that  satisfies  Eq.  (3.1.1)  with  c„  given.  From  Eq. 
(3.1.5): 


n 


Cn  e  2 


(3.4.6) 


STEP  6:  Compute  m  *(^)  and  the  wavelet  coefficients  dn 
From  Eq.  (3.3.2) 


mci^)  = 


^(^2;r/)| 

Vi!  |0(£t>f2;rZ)| 


(3.4.7) 


where  are  the  Fourier  coefficients.  They  can  be  obtained  via  a  Fourier  series  expansion; 


^  me  if) 


(3.4.8) 


STEP  7:  Construct  the  orthogonal  wavelet  y/*(.x) 
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From  Eq.  (3.3.8) 

-n-  1 


Remark:  Recall  from  Eq.  (3.3.6): 


i{x)  =  VIS  dnp  (2x-n) 


(3.4.10) 


Substituting  Eq.  (3.1.12)  into  Eq.  (3.4.10)  yields  an  alternative  way  to  construct  the  orthogonal 
scaling  function: 


0\x)  =  VIS  dn\^  am  (f>(2x-m-n)] 

n  m  -I 


(3.4.11) 


It  is  noted  that  the  construction  of  an  orthogonal  scaling  function  <(>  (x)  by  using  Eq.  (3.4.11) 
should  be  identical  to  that  of  using  Eq.  (3.4.5); 

STEP  8:  Fourier  Transform  of  <{>  (x)  and  \l/*{x) 

From  Eq.  (3. 1.9),  the  Fourier  transform  of  the  orthogonal  scaling  function  is  given  by: 


0  {CO)  = 


0(tU) 


|(;)(GH-2;r/)| 


=(X  ^ 


0(0)) 


(3.4.12) 


with  a.  defined  in  Eq.  (3.1.10b),  Similarly  for  Eq.  (3.2.18),  the  Fourier  transform  of  the 
ortho 2onal  wavelet  is: 


w\co)=^J^d.„.i  (-1)-"-!  0  if) 

V2  n 


(3.4.13) 


Example  I.  The  Piecemse  Linear  Spline  Function 

To  provide  a  simple  example,  we  choose  (t)  to  he  the  piecewise  linear  spline  function. 
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|l+x  -1<j:<0 

0(x)  =  \l-x  0<x.<l 
\0  otherwise 
as  plotted  in  Figure  3.1a. 

STEP  1:  This  function  satisfies  [see  Fig.  3.1a;  Eq.(3.1.1)] 

=  ^(Ix-l)  +  ^2x)  +  ^^2x+\y,  c.i  =  Co  =  Cl  = 

STEP  2:  From  Eq.  (3.4.3),  bn  are  computed  as  (two-point  Gaussian  integration  is  used). 

i>.,=r,6„  =  2.,i,,=i 

Thus  the  normalization  coefficient,  Eq.  (3.4.2a),  is  given  as: 

X  \i{Q}^27:lf  =  X  +  6 "  3  ^  6  “ 

;  n 

STEP  3:  The  orthogonal  expansion  coefficients  «„  as  defined  in  Eqs.  (3.4.4a,b)  are  tabulated  in 
Table  3.1.  / 

Table  3.1  Orthogonal  expansion  coefficients 
Note  that  a„  =  cc.„.  and  the  high  order  coefficients  cx„  are  neglected. 

STEP  4:  The  orthogonal  scaling  function  can  be  computed  via  Eq. (3.4.5).  It  is  shown  in  Figure 
3.1b. 

STEP  5:  The  computation  of  md^)-  From  Eq.  (3.4.6) 
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-f  =  X(. 

2  VI  2VI 


4re'f  +  4^  + =  i- (1  +  cos 


VI  2VI 


STEP  6:  The  compulation  of  m*(^).  From  Eq.  (3.4.7) 


me*  (y)  = 


V? 

'/l' 

^  2 

0(Q>f2;r/) 

=  1(1  +  cos  2 


./l  +  l 

(D)l  3  6 

^  /  /“TS  T" 

V 


cos 


1  +  1  cos  CO 
3  6 


-n.e  wavelet  coeff.cients  d,  can  be  computed  according  to  Eq.  (3.4.8)  and  they  are  given  in  Table 


n 

da - _____ - 

0 

0.8176464057010934 

1 

0.3972970881341911 

2 

-6.9 10098674164672E-02 

3 

-5.194534808183847E-02 

4 

1 .697 1 047 89387069E-02 

5 

9.9905954441 83770E-03 

6 

-3.88326225090555 lE-03 

7 

-2.20 195 1238397208E-03 

8 

9.2337 1005487 1564E-04 

9 

5.116360226930585E-04 

10 

-2.242963267262629E-04 

Table  3.2  Wavelet  coefficients,  d„  for  linear  spline  function 


Note  that  dr,  =  The  coefficients  d,,  are  nonzero  for  all  n  and  the  higher  order  coefficients  are 
neglected. 

STEP  7;  Using  expression  Eqs.  (3.4.11),  the  orthogonal  wavelet  is  shown  in  Figures  3.1c. 

STEP  8:  The  Fourier  transform  of  <p  (x)  and  y/  (x).  For  the  linear  spline  function  <t>{x),  the  Fourier 
transform  is  given  by  (see  Figure  3.  Id): 
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CO 

\  y  / 


Thus  the  Fourier  transforms  of  (p  {x)  and  i?*(x)  are  given  by  (see  Fig.3.1e  and  3. If,  respectively): 


(p  (0))  = 


sin^l 

2 

1  1  ■ 

1-^1 

2 

i  f  1 

1 

ta  j 
2  1 

't 

/l  +  lcos(0  „  1 

'  3  6 

VH 


COS 
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Example  2.  High  Order  Spline  Functions 

In  this  example,  we  present  the  various  coefficients  generated  for  the  cubic  spline  and  5th 
order  spline  function  (Tables  3.3  -  3.7).  The  spatial  and  Fourier  transform  space  for  the  scaling 

function  and  wavelet  are  depicted  in  Figures  3.2  and  3.3. 

Comparing  the  frequency  spectrum  of  the  linear,  cubic  and  5th  order  orthogonal  scaling 
functions  and  wavelets  (Figs.  3.1e  and  3.1f,  Figs.  3.1e  and  3.1f  and  Figs.  3.1e  and  3.1f. 
respectively),  it  can  be  seen  that  the  higher  order  scaling  function  approximates  a  rectangular  low- 
pass  filter,  whereas  the  higher  order  wavelets  have  better  high-pass  filter  properties.  It  is  noted 
that,  the  "side-lobes"  are  decreased  when  the  order  of  the  spline  increases. 


(bix)  (cubic  spline  function) 

(bix)  (5th  order  sphne  function) 

1(a:-h2)3,  -2<x<-1 

|-x2(l+|),  -1<X<0 

|-x2(l-|),  0<X<1 

4  (^  -  2)^  1  <  X  <  2 

6 

^(3  +  a:)5,  -3<a:<-2 

-J— (51  -75x-210x2-  150x3-45x^-5x5),  -2<x<-l 

120 

-L  (33  -  30x2+ 15x4  +  5x5),  -l<;c<0 

60 

(33 . 30x2  +  15x4  -  5x5),  o<x<l 

60 

-i-(51  +75x-  210x2  +  150x3-45x4  +  5x5),  1  <x<2 

120 

Table  3.3  The  cubic  and  5th  order  spline  functions 


The  coefficients  c„ : 


n 

Cn  (cubic  sphne  function) 

c„  (5th  order  sphne  function) 

-3 

0 

1 

_ 32^1 - 

1 

-2 

8V7 

_ \fdl - 

AL. 

-1 

_ 32dl _ 

A 

_3_ 

0 

4V7 

_ m. _ 

1 

15 

1 

9V7 

_ _ 22dl _ 

1 

3 

2  . 

8VT 

16VT 

3 

0 

1 

_ _ 32^1 _ 

Table  3.4  Coefficients  c„ 


The  coefficients  b„  : 


n 

bn  (cubic  spline  function) 

b„  (5th  order  spline  function) 

-4 

0 

0.393925565271235 

-3 

1.984126820776084E-04 

0.243960287323605 

.9 

2.380952406175427E-02 

5.520202023558386E-02 

-1 

0.2363095229476177 

3.8238786597 12875E-03 

0 

0.4793650806171007 

5. 1006093373373 19E-05 

1 

0.2363095229476177 

3.8238786597 12875E-03 

2 

2.380952406175427E-02 

5.520202023558386E-02 

3 

1.984126820776084E-04 

0.243960287323605 

4 

0 

0.393925565271235 

Table  3.5  Coefficients  b„ 


The  orthogonal  expansion  coefficients  a„  : 


n 

ccn  (cubic  spline  function) 

(5th  order  spline  function) 

0 

1.96976165984483 

3.21252767040872 

l' 

-.672430465527351 

-1.67129203556098 

2 

.2687042287778821 

.8693729316036274 

3 

-.118519934970734 

-.476253034889491 

4 

5.519145833683330E-02 

.2724049771107151 

5 

-2.652033542024175E-02 

-.160669074515415 

6 

1.2998 16380236 143E-02 

9.681 114156185701E-02 

7 

-6.457490594044601E-03 

-5.92 19828206672 19E-02 

8 

3.239862395307044E-03 

3.662439553866464E-02 

9 

-1.637775639515686E-03 

-2.2837729142741 13E-02 

10 

8.328359233600000E-04 

1.4332 19729959393E-02 

Table  3.6  Orthogonal  expansion  coefficients 


37 


The  wavelet  coefficients  d„  : 


n 

dn  (cubic  spline  function) 

dn  (5th  order  spline  function) 

0 

0.7661300483768689 

.7472336463549821 

1 

0.4339226328986687 

.442463060087634 

2 

-5.02017 19928895 17E-02 

-3.701993636800858E-02 

3 

-.11003701684689 

-.129268688614172 

4 

3.208089410823325E-02 

2.947409035770669E-02 

5 

4.206834992601 178E-02 

6.131245571925825E-02 

6 

- 1.7 1763 1339398893E-02 

-2.10061 1348868239E-02 

7 

-1.798231987132275E-02 

-3.252003030074334E-02 

8 

8.685293562579883E-03 

1.40095184725 1998E-02 

9 

8.201476481521336E-03  • 

1.820860315855871E-02 

10 

-4.353838725200284E-03 

-9.049 1 54645280076E-03 

Table  3.7  Wavelet  coefficients 


Figure  3.2  (a)  The  cubic  spline  function,  (b)  the  orthogonal  scaling  function,  (c)  the  wavelet,  (d) 
the  Fourier  transform  of  cubic  spline  function,  (e)  orthogonal  scaling  function,  (f)  orthogonal 

wavelet  (e) 


4.  REPRODUCING  KERNEL  METHODS  (RKM) 

4.1  Relations  between  Convolution  and  Reproducing  Kernel  Methods 

Given  a  sufficiently  smooth  function  u(x),  we  wish  to  construct  a  window  function 
<p(x)  e  which  is  also  called  a  reproducing  kernel,  such  that 


(p(x  -  y)  u(y)  dy 


(4.1.1) 


For  computational  efficiency,  (fix)  is  chosen  so  that  it  is  non-zero  only  in  a  compact  support,  i.e. 

for  every  B(y}, 


<pix  '-  y) 


l>0  xe  £(y) 
I  0  x€  B(y) 


Equation  (4.1.1)  is  can  also  be  expressed  as  a  convolution  integral,  ie,  (p(x)  *  u(x)  -  u^{x).  From 
the  convolution  theorem,  the  spatial  convolution  uHx)  is  equivalent  to  a  multiplication  in  the 

Fourier  domain. 


The  physical  meaning  of  the  convolution  Eq.  (4.1.1)  can  not  be  explained  clearly  either  verbaUy  or 
pictorially.  However,  it  becomes  more  transparent  when  we  study  its  Fourier  transform 

counterpart. 

In  order  to  control  the  frequency/wave  number  content  or  the  ‘scale’,  we  introduce  a 
scaling  parameter  (dilation  parameter  or  refinement  parameter  are  other  names)  a  >  0  so  that  Eq. 
(4.1.1)  becomes 


lA  ix)  =  PaUix)  s  U^(x  -  y)u{y)dy 


(4.1.3) 


where  the  projection  operator.  Pa,  is  defined  as: 
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PaU{x)  = 


oo 


{x  -  y)u{y)dy 


(4.1.4) 


and  a(y)^r*h(y)  with  r  =  2^  -2-1, 0,1,2,- -,00)  and /jfyj  is  the  nodal  spacing 

function.  Based  on  the  resolution  of  the  projection  operator,  ,  a  hierarchical  representation  ot 

the  function  u{x)  is  defined  as 

U{x)  =  lim  PyaU{x)  3  ...  PaU{x)  3  PlaUix)  3  />4.w(x)  3  ...  3  lim  =  {0}  (4.1.5) 


Now,  we  introduce  a  wavelet  projection  given  by 

Qia  u{x)  = 

where  i/iafx  -  >’)  is  the  2a  scale  wavelet.  It  is  defined  by 


y/iaix  -  y)  u(y)  dy 


(4.1.6) 


Vlaix  -  y)  =  (paix  -  y)  -  <ha{x  -  v) 


(4.1.7) 


Rewriting  Eq.  (4.1.4)  together  with  Eq.  (4.1.7)  yields 


PaU{x)  =  j  [(paix  -  y)  -  (piaix  '  .v)  +  02aix  '  .>’)]«(>•)  dy 

I  Viaix  -  vOwCv)  dy 

=  PzaUix)  +  QlaUix)  =  P^aHix)  +  04a«(^)  +  QlaUix)  =  ...  CtC.  (4. 1.8) 

Equation  (4.1.8)  illustrates  the  framework  of  multiresolution  analysis,  which  can  be  described  by 
sequences  of  nested  closed  subspaces  defined  in  Eq.  (4.1.5).  If  a  complementary  projection 
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operator,  is  defined  as  in  Eq.  (4.1.6),  then  the  higher  scale  projected  solution,  ,  can  be 
represented  by  the  sum  of  the  next  lower-level  scaling  projection  and  wavelet  projection.  The  Qia 
projection  can  be  viewed  as  a  “peeled  off’  scale  or  the  “rate  of  variation”  of  The  wavelet 
projection  Eq.  (4.1.6)  together  with  the  recursive  two-level  decomposition  Eq.  (4.1.8)  constitute 
the  backbone  of  the  wavelet  RKM. 

4.2  Reproducing  Kernel  Method  In  One  Dimension 

Reproducing  Condition 

Consider  the  reproducing  equation  as: 

u^-{x)=  {u{y)(l>a{x-y)dy  (4.2.1) 

J_oo 


The  Taylor  series  expansion  for  u(y)  about  x  is: 


u{y)  =  u(x)  -(x-  y)u  (x)  +  ~ 


(/i  +  l)! 


(4.2.2) 


nl 


where  .x  =  .x  +  d{y  -  x)  and  0<  6  <\.  Substitute  Eq.  (4.2.2)  into  Eq.  (4.2.1),  which  yields 


u^'^ix)  =  u{x)mQ{a,x)-u  (x)mi{a,x)  +—{x)m2ia,x)  + 


ft  ^  1 

n!  n\ 


(4.2.3) 


where  the  kth  moment  of  the  window  function  is  defined  as: 


m^.(rt..x)=  =  0’  ^ .  ” 


(4.2.4) 
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Apparently,  if  we  want  to  reproduce  u^^{x)  correctly  up  to  nth  order,  the  following  reproducing 
conditions  need  to  be  satisfied 

m„(a,x)  =  5oa  a  =  0,  1,  2, n  (4.,u.5) 

An  arbitrary  choice  of  a  window  function,  however,  can  not  guarantee  the  satisfaction  of  all  the 
reproducing  conditions  Eq.  (4.2.5).  Therefore,  we  define  a  modified  window  function  so  that  all 
the  reproducing  conditions  can  be  satisfied.  In  general,  this  modified  window  function  can  be 

expressed  as: 

n 

'^^{x\x-y)  =  Y,^kM{x-yf<Pa[x-y)  =  P'^{x-y)b{a,x)(p^{x-y) 


=  Ca{x\x-y)(pa{x-y) 


where 


P^{x-y)  =  \l{x-y) . (x-y)" 


(4.2.6) 


(4.2.6a) 


and 


b'^{a,x)  =  lbo{a,x),b^{a,x),...,b„{a,x)]  (4.2.6b) 

Note  that  a  polynomial  type  vector  is  adopted  to  the  basis  function  in  Eq.  (4.2.6a).  However, 
it  is  not  the  only  option  for  the  basis  function  in  RKM.  In  general,  the  basis  function  can  be 

extended  to  any  independent  functions.  The  product  of  P  (x  — _v)  and  b{a,x),  Cq(x,x  y),  is 

called  the  correction  function.  The  coefficients  bi,{a,x)  are  solved  from  the  reproducing 
conditions.  By  this  definition,  the  moments  of  the  modified  window  function  are  written  as 


+  00 

/■«;.(«. x)  =  |(x  -  yf  0a(x;x  -  y)<iy 


k  =  0.  1,  2...,  n 
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=  f(jc-vf  '^bj{a,x){x  -  yY  (t>a{x  -  y)  dy 

L>=o 

=  bo{a,x)mi,{a,x)  +  b,{a,x)m,^i{a,x)+...+b„^^^^ 


or  in  matrix  form 

m{a,x)  =  M{a,x)b{a,x) 

where  M(a,  x)  is  the  moment  matrix: 


(4.2.8) 


M{a,x) 


mo(a,x) 

mfa,x) 

mfa,x) 

m2{a,x) 

■  "i„(a,x) 

m„{a,x) 

m„+i(a,x) 

•  m2„(fl,x) 

(4.2.9) 


To  solve  the  coefficients  we  apply  the  reproducing  conditions  Eq.  (4.2.5)  to  the  modified 

window  function.  That  is. 


[/7iQ(a.x),mi(a,x),---,m„(a,x)]  -[l,0,---,0]  P  (0) 

Equations  (4.2.8)  and  (4.2.10)  give: 

M{a.x)b{a,x)  =  P{0) 

Then,  the  coefficient  vector  b{a,x)  can  be  obtained  by: 

bia.x)  =  M-Ha,x)PiO) 

Reproducing  Condirions  for  the  First  Derivative 

The  first  derivative  of  the  reproducing  equation  Eq.  (4.2.1)  is  defined  as: 


(4.2.10) 


(4.2.11) 


(4.2.12) 


—u^^{x)=  f  i<{v)4-^a^^-y'>^y 
(lx  J-~  ■  dx 


=  {i4{y)D^<l)^{x-y)dy 

J  — oo 


(4.2.13) 
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Substituting  Eq.  (4.2.2)  into  Eq.  (4.2.13)  yields 


—u^‘‘{x)  =  u{x)mQi{a,x)-u  {x)m^_i{a,x) 
dx 


2!  n\ 


(n  +  1)! 


(4.2.14) 


where  is  defined  as: 


m„.i(a,x)=  \{x-yfD\{x-y)dy 
J—oo 


a  =  0,  1,  2, ...,  n 


(4.2.15) 


By  examining  Eq.  (4.2.14),  the  reproducing  conditions  up  to  nth  order  for  the  first  derivative  are: 

=  a  =  0,  1, 2 . n  (4.2.16) 


nir 


In  oeneral,  if  we  define  the  7th  derivative  of  the  reproducing  equation  as: 


u{y)D^<t)aix-y)dy 


(4.2.17) 


where  ‘D’  denotes  a  differential  operator,  then,  the  reproducing  conditions  for  the  7th  derivative 
can  also  be  obtained: 


mf^y{a,x)  =  (-l)“a!5„y 


(4.2.18) 


where 


y{a,x)=  f  (-1^  ~  .v)°'  D^0a^x  -  y)dy 

J  —oo 


(4.2.19) 
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Example  1:  From  Eq.  (4.2.6),  to  reproduce  the  first  derivative  with  up  to  second  order  accuracy, 
the  modified  window  function  is  obtained  in  the  form  of: 


0^{x\x-y)  =  \bo{a,x)  +  bi{a,x){x-y)  +  b2{a,x){x-yf'\(pa{x  .v) 


(4.2.20) 


Applying  the  modified  window  function  to  the  reproducing  condition  for  the  denvative  and 
by  introducing  the  reproduction  condition  ni(^{a,x)  =  yields: 


/Tlo 

rrii 

m2 

mo 

mi 

m2 

mi 

m2 

m3 

mi 

m2 

m3 

mj 

m3 

m4 

m2 

m3 

m4 

h 

bi 

bo 

bi 

Pi 


(4.2.21) 


Combining  Eq.  (4.2.11)  and  Eq.  (4.2.21)  gives 


*  M{a,x) 

0 

b{a,x) 

'p(oy 

M'{a,x) 

M{a,x)_ 

b  (a,x)_ 

0 

(4.2.22) 


Taking  the  first  derivative  on  both  sides  of  Eq.  (4.2. 11),  yields. 
M  {a,x)b{a,  x)  +  M{a,  x)b  {a,  x)  =  0 


(4.2.23) 


This  shows  that  the  reproducing  conditions  for  the  first  derivative  can  also  be  achieved  by  taking 
the  derivative  of  the  reproducing  conditions  directly.  That  is: 


ma.Mx)  =  D\M{a,x)b{a,x)\  = 

Recalling  the  reproducing  conditions  for  the  7 th  derivative  is  given  as: 

Tna,y{o^x)  =  {-\)^  a\  5  ay 


(4.2.24) 


(4.2.25) 


li  can  be  shown  that  these  conditions  simply  imply: 
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D^m(a,x)  =  D’'[M(a,x)b(a,  a:)]  =  0 


(4.2.26) 


4.3  Orthogonal  Conditions  for  RKM 

We  define  a  new  correction  function  associated  with  the  reproducing  and  orthogonality 
conditions  for  a  discrete  system  with  uniform  spacing,  as: 


Q{x,xj)  =  + 


X-Xj 


,X-Xj. 


(4.3.1) 


where  ‘n’  is  the  order  of  accuracy  of  the  reproducing  conditions  and  ‘m’  is  the  number  of  nodes 
covered  by  the  support  of  the  window  function  ^(-) .  Hence,  a  new  set  of  coefficients  Pa  ’ 
to  be  determined  from  both  the  reproducing  and  orthogonality  conditions. 

/ 

Reproducing  Conditions 

From  Eq.  (4.2.5),  for  0*(^~  — )  to  reproduce  the  solution  up  to  nth  order  requires: 

n+m 

J  =  . n 

i=0 

There  are  n+1  equations  with  n  +  m  +  l  unknowns.  In  order  to  obtain  a  complete  set  of  /3„  ( a  = 
0.  1 . n+1 .  n+m  ),  another  m  equations  can  be  acquired  from  the  orthogonality 

conditions. 


Orthogonalin-  Conditions  for  Scaling  Function 

Instead  of  solving  for  the  filter  coefficients  a„,  the  same  orthogonality  conditions  can  be 

applied  to  solve  the  correction  coefficients  j3„.  Furthermore,  with  the  property  ot  compact  suppon 

for  the  window  function,  only  the  nodes  covered  by  the  support  need  to  be  considered  in  the 

oithogonality  condition.  That  is: 
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oo 

J  (l>*  {x)0*  (x  -  k)  =  5qi^ 


(m-1) 


(4.3.4) 


Substituting  Eq.  (4.3.2)  into  Eq.  (4.3.4),  yields: 


~  n+m  n+m  ,  x/  /■  \  \ 


i=0  ;=0 


(4.3.5) 


where  k 


(m-1)  _i  0  1  Note  that  Eq.  (4.3.5)  implies  m  equations.  With  a 

/C  —  ?  •  •  •  ^  1,  u,  i,  . . . ,  •  ^ 

2  ^ 

total  of  n+m+7  equations,  an  orthogonal  window  function  can  easily  be  obtained  by  solving 
n+m+1  coefficients  /3„. 

Next,  we  will  introduce  a  new  set  of  correction  coefficients,  which  enforce  the  original 
window  function  to  satisfy  the  wavelet  orthogonality  conditions.  An  orthogonal  wavelet  function 
with  correction  coefficients  is  defined  as; 


*,X-Xj^ 

¥  (: - )  = 


To  +  yi(^^)+ 


a 


X  -_Xj 
a 


(4.3.6) 

a 


where  k 


support 


(^-^)  -10  1  and  ‘m’  is  the  number  of  nodes  covered  by  the 

2  .  2 

of  the  window  function  Again,  these  coefficients  r„  can  be  solved  from 


orthogonality  conditions  for  the  wavelet. 

Orthogonality  Condition  for  Wavelets 

The  requirement  for  an  orthogonality  between  wavelet  functions  and  scaling  functions  is; 


J  ¥*{x)(l)*{x-k)  =  0  k---  -  .  1,0,1 . 


(m-1) 


Subsiiiutina  the  definitions  of  V^*(x)  and  0  (x-k)  into  Eq.  (4.3.7),  yields. 


(4.3.7) 
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~  m  n+m  .  ,  \j 

IS; 

1=0  j=0 


dx  =  0 


(4.3.8) 


Where  k  =  -^'^  —  ..-1,0,1,.../^  Note  that  Eq.  (4.3.8)  also  implies  m  equations. 

2  2 

Finally,  to  accomplish  the  orthogonality  of  the  window  function  and  wavelet  function,  the 
determination  of  the  coefficients  necessary.  Although  the  coefficients  Ya  are 

entangled  with  in  Eq.  (4.3.8),  they  can  be  treated  as  a  set  of  linear  equations  once.  is 
obtained.  However,  it  is  still  difficult  to  obtain  closed-form  solutions  for  /3„  with  a  set  of  non¬ 
linear  equations,  Eq.  (4.3.5).  A  numerical  solution  will  be  the  possible  choice  for  obtaining  these 
coefficients. 


4.4  Two  Dimensional  (2D)  RKM 

By  a  direct  extension  from  one  dimension,  a  general  2D  window  function  is  defined  as; 


A  special  2D  window  function  can  be  represented  by  the  product  of  ID  window  functions.  Thus, 

The  basis  function  for  2D  is  given  in  the  form  of: 

P^{x)  =  [1,  P( {xl  pI (X), ...,  pI (x)]  (4-4. 


where 


pI (X)  =  [xf ,  .xrf ■'.X2 . ,  x\ ]  kth  order  polynomials 


(4.4.4) 


Reproducing  Conditions  in  2D 

The  reproducing  equation  for  2D  is  now  defined  as: 
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m&o 

u^‘(x)=  \u{y)0a(x-y^^y 

J  — oo 


(4.4.5) 


Substituting  a  2D  Taylor  series  expansion  into  Eq.  (4.4.5)  yields; 

u^-ix)  =  m^(a,x)uix)-^mio{a,x)Ux^ix)  +  mQiia,x)u^^ix)^ 


+  i  m2o(a,^)^  +  2mn(a,x)^j^  +  mo2(«.^)^^2j“(^)+^ 


(4.4.6) 


where  a  2D  moment  equation  is  defined  as: 


(a,Xi,X2)=  J  (Xi  -  yiY  {X2  -  y2y  ~  y^‘^y 


0<l,j<k 


(4.4.7) 


From  Eq.  (44.6),  it  can  be  seen  that  to  reproduce  the  solution  correctly  up  to  «th  order  m  2D,  the 
reproducing  conditions  require  that; 


,{a,x)  =  5qj5q 


0<  I,  J<n 


(4.4.8) 


To  satisfy  the  reproducing  conditions,  the  modified  window  function  for  2D  is  assumed  to  be: 


0„{x)  =  [p'‘' {x  - y)b{a,x)]  <^„(x) 


(4.4.9) 


where 


b(a,xf  =  [^7oo("»^)'  ...,b„Q{a,x) . bQ„{a,x)] 


i4.4.10) 


Then,  the  moments  of  the  modified  window  function  in  2D. are  denoted  as; 
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m,j{a,x)=  r  {x,-y,y{x2-y2)\P^{x-y)b{a,x)]<t>aix-y)dy 


oo  n 

=  J  {xi-yO^ixi-yi^  ^ 


{x^-yxY{x2-y2Yhrs{a,x) 


n 

'^mp^{a,x)b„{a,x)  p  =  r  +  Imdq  =  s  +  J 

r,s=0 


or  in  matrix  notation; 


m{a,x)  =  M{a,x)b{a,x) 
where  M{a,x)  is  the  moment  matrix  in  2D,  thus: 


M{a,x) 


moo 

mxo 

"Joi 

•  m„o  ■ 

■■  mon 

m\o 

m2o 

mil 

mi„ 

mo\ 

mil 

mo2 

mon+l 

m„Q  : 

tnor,  mxn  rn^n+\  2n 


Enforcing  the  reproducing  conditions  Eq.  (4.4.8)  on  Eq.  (4.4.11)  yields: 


M{a,x)b{a,x)  =  PiO) 

Then,  the  coefficient  vector  b{a,x)  can  be  computed  by  matrix  inversion 
b{a,x)  =  M~Ha,x)PiO) 


(4.4.11) 


(4.4.12) 


(4.4.13) 


(4.4.14) 


(4.4.15) 
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4.5  Multi-Dimensional  RKM 

To  present  RKM  in  multiple  dimensions,  multi-index  notation  is  adopted. 

Multi-index  Notation 

We  will  use  multi-index  notation  in  the  following,  a  (a  =  [aj,a2’' '  ’  ^  multi-index 

if  «,  (/  =  1,2,  -  •  W)  are  non-negative  integers  with  the  foUowing  properties  for  this  multi-index  a  : 


N 


1=1 


a\  =  ai\a2'-cci^\ 


(4.5.1a) 


(4.5.1b) 


We  then  define  two  symbols 


m 


1  ^2'  '"^N 


(4.5.1c) 

(4.5.1d) 


With  this  multi-index  notation,  the  Taylor  series  expansion  in  terms  of  degree  n  followed  by  a 
remainder  can  be  written  as: 


|ai=0  ^  ’  |a|=n+l 


(4.5.2) 


To  understand  how  this  multi-index  notation  works  for  the  Taylor  series  expansion,  a  3D  example 
is  shown  next. 

Example  1.  3D  Taylor's  formula  using  the  multi-index  notation. 

In  3D.  the  properties  for  multi-index  notation  are  (  Eq.(4.5.1) ): 


i  \a\  a  ■«,  -I- «;  -I-  a? 


a!  =  a, !  a-. !  a?  1 


(4.5.3a) 

(4.5.3b) 
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iii  {x-yf  {4.53c) 

iv  D‘^u{x)  =  dy^-d^;u{x)  (4.5.3d) 

Table  5.1  lists  all  the  parameter  values  for  the  multi-index  and  shows  how  this  multi-index  notation 
works  for  3D  Taylor’s  formula. 


1«1 

«i 

^2 

«3 

a! 

D“u(x) 

^  Ul-£)“«(A:)(x-y)“ 
a\ 

l«l=o 

0 

0 

0 

0 

1 

1 

u(x) 

uix) 

1 

1 

0 

0 

1 

u(x) 

■*1 

-(ATi  -yi)^;,^U(X) 

0 

1 

0 

1 

(.>;2 

d^.u{x)  . 

-(j:2-y2  )d^u(x) 

0 

0 

1 

1 

(j;3-.V3) 

d,u(x) 

-{Xi-yi)d^u(x) 

2 

2 

0 

0 

2 

(?*U(X) 

i[(x,-yir<?x. 

0 

2 

0 

2 

(^2.- >”2)* 

+(  X2  ”  VS  )*  ^x-, 

0 

0 

2 

2 

(j;3-.V3  )* 

d;u(x) 

+(*3-y3)'^x, 

1 

1 

0 

1 

(■>;i->’l)(^2-.'’2) 

d^^d^u(x) 

+2(x,-y,  )(X2-y2)d^^d^ 

0 

1 

1 

1 

(Jt2-)i2)(A:3-y3) 

d^,d^u(x) 

+2(x2-y2){X2-y})dx,dx, 

1 

0 

1 

1 

(Xj-Vi  )(A:3-y3) 

d,d^u(x) 

+2(x,  -y,  )(.t3  -  y3  ^ 

3 

3 

0 

i 

6 

6 

dlu(x) 

0 

3 

0 

6 

(^2  - 

.dl^uix) 

+(•^2  '‘>’2 

0 

0  ' 

3 

6 

(J:3-y3)^ 

dlu(x) 

2 

1 

0 

2 

(Xj-Vi  )’(.V2-y2) 

+3(A:i-yi  riX2-y2)^l,^x, 

2 

0 

1 

2 

(.tl-yi)*(.t3-.V3) 

d;d^u(x) 

1 

2 

0 

2 

(j:,  -Vi  )(Ar2  -  .V2  r 

+3(  —  >’1 )( AS  —  vs  )  3^2 

0 

2 

1 

2 

(^2-.'’:)‘(-'f3">'3) 

did,u{%) 

4‘3(  X2  ■“  vs  )"■  {  A3  —  >'3  )3x2  ^Xj 

0 

1 

2 

2 

(^2--''2)(-’^3“>’3)‘ 

d,,dlu(x) 

+3(  AS  ys )( A3  •”  y3  )‘^  ^X2 

1 

0 

2 

2 

(.V, -y,  K-V3-.V3)‘ 

+3(Ai-yi  )(A3-y3 

1 

1 

1 

1 

(ATi-y,  )(a:2-.V2  )(Jt3-.V3) 

d,d^d,u(x) 

+6(A^  ‘“  Vj  )(A2  “  3*2  ){A3  ^X’.  ^.tj  ) 

Table  5.1  Terms  in  3D  Taylors  Formula  in  multi-index  notation 


Miihi-Dimensional  RKM 

The  general  window  function  for  N  dimensions  is  defined  as: 
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Assume  an  A- dimensional,  n-order,  /-components  basis  function  is  given  as: 

P^(x)  =  [l,  P[ (xX  pI (X), ...,  pI (X)] 

where 

fJ(x)  =[xi",  Xi"~^X2, ...,  X^v-l^r'’  Av 

R€pwducing  Conditions  in  N-Dimcnsions 

The  reproducing  equation  for  multiple  dimensions  are: 

J_oo 

Substitution  of  Eq.  (4.5.2)  into  Eq.  (4.5.7)  by  using  the  modified  window  function 
gives: 

„^<!(x)=^  (9“M(x)m„(a,  x)+  ^  d^u{x)ma{a,  x) 

|al=0  |al=n+l 


(-1) 

where  Cf  is  a  multi- index,  d“u{x)  =  —D°u{x)  and 

GC . 

m«(a,x)=f  {X  -  -  y)dy 

J_e>o 

The  a  moment  of  the  modified  window  function  0«(x;x  -y)  satisfies: 


niaia,  x)  =  5^a 


(4.5.4) 

(4.5.5) 

(4.5.6) 

(4.5.7) 

’fl(x;x-y) 

(4.5.8) 

(4.5.9) 

(4.5.10) 


From  previous  analysis,  the  reproducing  conditions  in  Eq.  (4.5.10)  imply 


M{a,  x)b{a,  x)  =  PiO) 


(4.5.11) 


wheri  M{a,  x)  consists  of  moments  of  the  window  function  of  0a{x;x-y),  that  is: 


^00 

o 

o 

••  M„o' 

Moi 

Mil 

Mnl 

M{a,x)  = 

Mq2 

Mon 

Mi„  -  • 

••  A/„„, 

(4.5.12) 


and  the  components  of  the  moment  matrix  are; 

M„^(fl,x)  =  j  {X  - yf  <Paix-,x  - y){x  -  yf  dy  0<|a|,l^l< 


(4.5.13) 


where  both  a  and  (3  are  multi-indices. 

Reproducing  Conditions  of  Derivatives  in  N -Dimensions 

Extending  the  ID  analysis,  the  reproducing  conditions  of  the  yth  derivative  for  N 

dimensions  ^e  given  by: 

These  conditions  imply: 


D^[M{a,  x)b{a,  x)]  =  D'^{P{0)] 


or  in  the  matrix  form: 


(4.5.15a) 


M  0  0 

DM  M 

D-M  2DM  M 


0 

0 

0 


Db 

D^b 


'p(oy 

0 

0 


(4.5.15b) 
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where  Cy 


y- 

(y-m)!m! 


{m  =  0, 1, 2,  •••.  y) 


(4.5.16) 


5.  REPRODUCING  KERNEL  PARTICLE  METHODS  (RKPM) 

The  reproducing  kernel  particle  method  (RKPM)  is  the  discrete  form  of  RKM.  For  the 
discrete  system,  the  reproducing  conditions  can  be  obtained  following  an  extension  ot  the 
procedure  developed  in  the  continuous  case. 

5.1  Reproducing  Kernel  Particle  Method  In  One  Dimension 

Reproducing  Condition 

The  reproducing  equation  in  the  discrete  system  is  given  as: 

np 

lAix)  =  '^u{Xj)<l>a(x-Xj)AXj 

where  np  is  the  number  of  the  particles  in  the  discrete  system.  We,  then,  define  the  numencal 
moments  o/  the  window  function  as: 

np 

m^{a,x)  =  '^(x-xjf  (t>a(x-Xj)Axj  /:  =  0,  1,  n  (5-1-2) 

;=i 

Recall  from  Eq.  (4.2.3),  to  reproduce  (x)  exactly  up  to  nth  order,  we  require. 

ma{a,x)  =  doa  a  =  0,  1, 2, ....  n  (5.1.3) 

To  satisfy  the  reproducing  conditions  Eq.  (5.1.3),  a  modified  window  function  in  the  discrete 
system  is  assumed  to  be  in  the  form  of: 

d„(x:x  -  Xj)  =  '^bk{a,x){x  -  Xjf  ((>a(x  -  Xj)  =  P^(x  -  xj)b{a,x)(p^(x  -  .x^) 
k=0 
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=  Ca(x-,x-Xj)(t)^(x-Xj)  (5.1.4) 

The  basis  function  in  Eq.  (5.1.4)  is  chosen  as  a  polynomial  type  vector.  The  coefficients 
bk{a,x)  can  be  solved  by  applying  the  modified  window  function  to  the  reproducing  condition. 

That  is: 

M{a,x)b{a,x)  =  P{0)  (5.1.5) 

where  M(a,x)  is  the  numerical  moment  matrix: 


’/no(fl,x) 

mi(a,x) 

•••  m„{a,x) 

mi(a,x) 

/^(a,x) 

_m„{a,x) 

m„+i(a,x) 

^2ni^'X)_ 

np 

- 

7=1 

^j)^a(x-Xj 

! 

T 

Then,  the  coefficient  vector  b{a,x)  can  be  obtained: 
bia,x)  =  M~^{a,x)P{0) 


(5.1.6) 


(5.1.7) 


This  RKPM  can  be  viewed  as  a  modified  convolution  formula  for  a  finite  domain  and 
nonuniformly  sampled  system.  The  modified  window  function,  is  the  low-pass 

filter  for  this  finite  domain,  which  can  be  viewed  as  a  nonuniformly  sampled  system.  In 
computational  mechanics,  we  use  this  to  define  the  shape  function,  such  that 

Nj{x)  =  Ca{x,Xj)0a(x-Xj)^x  shape  function  (5.1-8) 

np 

and  n^jx)  =  ^  u{Xj)Nj{x) 


reconstruction 


(5.1.9) 
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Reproducing  Conditions  for  the  Derivative 

The  7th  derivative  of  the  discrete  reproducing  equation  is  defined  as; 


=  '^u{Xj)D^(t>a(x-Xj)AXj 


(5.1.10) 


where  '£>’  denotes  a  differential  operator,  then,  the  reproducing  conditions  for  the  7th  derivative 
can  also  be  obtained; 


(5.1.11) 


where 


=  ^{x-xjf  D^tp^ix  -  Xj)Axj 


(5.1.12) 


Example  1:  From  Eq.  (5.1.4),  to  reproduce  the  first  derivative  up  to  the  first  order  accuracy,  the 
modified  window  function  is  assumed  to  be  in  the  form  of; 


(p^{^x;x  -  Xj^  =  ^boia^x)  +  bi{a,x){^x  -  Xj^y)a[x  Xj^ 


(5.1.13) 


Enforcing  the  reproducing  conditions  /n„(a,  j:)  =  <5oa  and  the  reproducing  condition  for  the 
derivative  m.Qc  -[{a,x)  =  —5a\  on  the  modified  window  function  yields. 


/riQ  mj  mQ  m^  tn 

my  m2  my  m2\ 


(5.1.14) 


or  in  the  matrix  form 


M  {a.x)b{a,x)  +  M{a,x)b  {a,x)  —  0 


(5.1.15) 
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From  Eq.  (5.1.15),  it  can  be  seen  that  the  reproducing  conditions  for  the  first  derivative  can  also  be 
obtained  direcdy  by  taking  the  derivative  of  the  reproducing  conditions.  That  is; 

D^[M{a,x)b{a,x)\  =  D^P{0)  (5.1.16) 

5.4  Two  Dimensional  (2D)  RKPM 

The  RKPM  can  be  easily  extended  to  multiple  dimensions.  A  general  2D  window  function 
is  assumed  of  the  form: 

^  .  X  ^  ^  (5-2.1) 


The  basis  function  for  2D  is  given  in  the  form  of; 
P^{x)  =  [l,  pI ix),  pI (X) . pI (X)] 


where 


kth  order  polynomials 


(5.2.3) 


Reproducing  Conditions  in  2D 

The  reproducing  equation  for  2D  is  now  defined  as. 


np 

U^‘{X)  =  (X  -  Xj  )AX  j 

;=i 

where  Axy  denotes  the  jth  nodal  area.  The  moment  in  2D  is  defined  as: 


(5.2.4) 


np 

m,j(a.x.  y)  =  ^  (x-x^)^y-yp^<^>„(x-Xy)Ax^  0<LJ<k 
;=i 


(5.2.5) 


where  k  =  0.  1 . n.  From  Eq.  (4.4.6),  it  can  be  shown  that  to  reproduce  the  solution  correctly 

up  to  /nh  order  in  2D,  the  reproducing  conditions  are  in  the  form: 
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m[j(a,x)  =  5qi5qj  0<I,J^n 

The  modified  window  fimction  for  2D  is  assumed  to  be; 


0^ {X-X  -Xj)  =  [p'^(x- Xj )b{a, X)]  0„ (x;x  - Xj) 


(5.2.6) 


(5.2.7) 


where 


bia,xf  =[bQo(a,x),bio{a,x),bQiia,x) . b„Q(;a,x) . bo„ia,x)] 


(5.2.8) 


Enforcing  the  reproducing  conditions  Eq.  (5.2.6)  on  the  modified  window  function  yields: 


M{a,x)b{a,x)  =  F(0) 


where  M{a,x) ,  the  numerical  moment  matrix  in  2D,  is  given  as: 


M{d,x) 


1^10 

ibol  ■■ 

•  ^nO  ■■ 

■■ 

^10 

dl20 

mil 

ibin 

dtoi 

mil 

mo2 

dtOn+l 

difio 

dlQn 

dlln 

dlOn+\ 

dh2n 

Then,  the  coefficient  vector  b{a,x)  can  be  computed  by  matrix  inversion 
b{a,x)  =  M~^{a,x)P{Q) 


(5.2.9) 


(5.2.10) 


(5.2.11) 


Reproducing  Conditions  for  Derivative  in  2D 

Extending  the  ID  analysis,  the  reproducing  conditions  of  7th  denvative  tor  2D  are 


np 

I 


(.x-.x, 


(  V  -  y^)“=  d'^'  d'^-  0a{x  -Xj)lSXj  =  (-1)“'  “-  «! !  02  ! 


(5.2.12) 
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Following  the  same  procedures  in  Examples  5.1,  it  can  be  shown  that  the  reproducing  conditions 
in  Eq.  (5.2.12)  simply  imply: 

D^m{a,  x)  =  x)b{a,  x)]  (5.2.13) 

5.5  Three-Dimensional  (3D)  RKPM 

The  general  window  function  for  three  dimensions  is  defined  as: 


A  basis  function  is  given  in  the  form: 

P^{x)=[\,P({x\P2(x\...,Pji{x)\ 

where  Plix)  are  complete  polynomials  of  order  k.  For  example,  P^  (x)  =  [x,y,..], 
p|'(x)  =  [.x^,y^,z^,xy,xz,yz],  etc. 

/ 

Reproducing  Conditions  for  3D 

Recall  from  2D  case,  the  reproducing  equation  for  3D  is  given  as: 


"P 

u^‘{x)  =  ^u{Xj)0a{x-Xj)Ax  J 

;=i 

where  Axj  denotes  the  jth  nodal  volume.  The  moment  in  3D  is  defined  as: 

np 

y. :)  =  y  (x-Xj)\y-yjy(z-z,f*Ax 
Recall  Eq.  (4.5.10),  the  reproducing  conditions  are  in  the  form: 


(5.3.3) 


(5.3.4) 


-  ^OI^OJ^QK 


(5.3.5) 
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From  previous  analysis,  the  reproducing  conditions  in  Eq.  (5.3.5)  imply 
M{a,  x)b{a,  x)  =  P(0) 

where  the  moment  matrix  M{o,x)  for  3D  is  given  by. 


M{a,x)  = 


^00 

M,o 

Moi 

Mil 

Mq2 

Mon 

M,n 

and  the  components  of  the  numerical  moment  matrix  are: 

np 

=  ^  (x-Xj)“0Jx-.x-Xj)(x-XjfAxj  0<|a|,|^|S* 


where  both  a  and  (5  are  multi-indices. 


(5.3.6) 


(5.3.7) 


(5.3.8) 


Reproducing  Conditions  of  Derivative  for  3  D 

From  2D  analysis  and  with  the  help  of  multi-index  notation,  the  reproducing  conditions  of 

the  }th  derivative  for  3D  are  given  by: 


ma.y{a,x)  =  {-\fa\6ya 

These  conditions  imply: 

x)b{a,  x)]  =  D^[PiO)] 


(5.3.9) 


(5.3.10) 


or  in  the  matrix  form: 


M 

0 

0 

0 

■  b 

'P(O) 

DM 

M 

0 

Db 

0 

D-M 

2  DM 

M 

0 

D^b 

— 

0 

CyD'^M 

c\d^'^m 

■■  C]M 

p'^b 

0 

where  C"' 


Y' 

(y-m)!m! 


(m  =  0, 1,  2,  •  •  •,  y) 


The  3D  computational  aspect  are  given  in  Appendix  A. 

6.  REPRODUCING  KERNEL  PARTICLE  METHODS 
IN  FOURIER  SPACE 


6.1  Preliminary 

From  the  Fourier  analysis,  the  convolutions  of  polynomials  ot  order  zero  (1), 
two  (x^)  with  a  window  function  can  be  expressed  as 


|l'(/)a(x-y)^iy  =  ^(0). 


j  y  -(paix  -  y)dy  =  4(0)  -  (O)  =  -  i^;(o) 


-1-00 


J  y^  -  y)dy  =  ^^^(0)  -  i2ax0'(O)  -  a4' 


(0) 


=  x4a(0)-i24;(0)-C(0) 
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(5.3.1  n 


(5.3.12) 


one  (x)  and 


(6.1.1) 


(6.1.2) 


(6.1.3) 


where  is  the  nth  derivative  of  evaluated  at  x=0 

By  induction,  the  convolution  of  a  polynomial  of  order  n  with  a  window  function  is  given  as 
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Tv"  .a(.v  -  v)A'  =  y  (-/)‘c;."-‘V4“'(o)  = 

J  S  ^--0. 


(6.1.4) 


and 


G"  = 


n; 


(6.1.5) 


6.2  Reproducing  Conditions 


In  order  to  reproduce  the  polynomials  correctly,  several  requirements  of  the  window 
function  need  to  be  established  in  the  Fourier  transform  domain  (^: 

1.  To  reproduce  [1],  we  need  ^(0)  =  1.  ^  (6.2.6a) 

2.  To  reproduce  [1,  a:],  we  need  0(0)  =  1  and  ^'(0)  =  0 .  ^  (6.2.6b) 

3.  To  reproduce  [1,  x,  x^],  we  need  0(0)  =  1,  0'(O)  =  0  and  0"(O)  =  0.  (6.2.6c) 

Generalizing  these  reproducing  conditions  to  nth  order  polynomials  yields  . 

/ 

To  reproduce  [1,  ;c . we  need  0(0)  =  1  and  =  0  for  A:  €  [l,n].  (6.2.7) 


These  reproducing  conditions  set  a  constraint  on  the  selection  of  a  window  function.  The  higher 
the  order  of  polynomial  requested,  the  flatter  the  window  function  must  be  at  <^=0  in  the  Founer 
transform  domain.  Accordingly,  an  ideal  low-pass-filter  window  function  (a  window  function 
with  all  the  derivatives  equal  to  zero  at  (^=0)  would  reproduce  a  polynomial  of  arbitrary  order. 


6.3  Reproducing  Conditions  in  Terms  of  Window  Function  Moments 

In  our  research,  we  are  not  interested  in  selecting  a  window  which  reproduces  a  given 
function,  but  in  designing  the  window  to  have  certain  filtering  properties.  Therefore,  an  equivalent 
expression  of  the  reproducing  conditions  is  needed  in  the  function  domain  (x).  First  of  all,  let  us 
define  the  moments  of  the  window  function  as 
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-^oo 

mk{a,x)=  j{x-  yf  0a{x-  y)dy  kth 


moment  of  window  function 


(6.3.1 » 


The  convolutions  defined  by  moments  of  a  window  function  are  given  as 


-1-00 

^{■0a{x-y)dy  =  mo{a.x) 


(6.3.2) 


-1-00 

^ y  '0a{^  ~  y)‘^y  ~  J[^  ~  6  ~ '^ai^  ~  y)^y  ~  xmQ{a,x)  -  mi{a,x) 


(6.3.3) 


+0O 

I  y^  “  y)^y  =  J  "  ^^6  ■  y) + 6  ~  >’)^]  ~ 


=  x^mQ{a,x)-2xmi{a,x)  +  m2{a,x) 

To  summarize,  the  convolution  of  an  nth  order  polynomial  is  of  the  form 


(6.3.4) 


J .v”  -  y)dy  = 

„  k=0 


(6.3.5) 


From  Eq.(6.3.5),  reproducing  conditions  similar  to  Eqs.  (6.2.7)  can  be  derived 

1.  To  reproduce  [1]  we  need  mo (fl,x)  =  1.  (6.3.6a) 

2.  To  reproduce  [1,  ;c]  we  need  mo(a,x)  =  1  and  m^{a,x)  =  0.  (6.3.6b) 

3.  To  reproduce  [1,  x,  x^  we  need  m^iyiyx)  =  1,  mi{a,x)  =  0  and  ini{a,x)  =  0  (6.3.6c) 

In  general,  the  reproducing  conditions  for  the  nth  order  polynomial  in  the  form  of  moments  are 
given  as 

To  reproduce  [1,  x . x"]  we  need  mo(a,:r)  =  1  and  m^.{a,x)  =  0  for  ^  e  [l,n]  (6.3.7) 
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Eq.(6.2.7)  and  Eq.{6.3.7)  can  be  related  by  the  moment  equations 
m,(.r)  =  /V0''''(O)  =  /'-0f  (0)  for  A:€[0,n] 


where  /  =  v-1  • 

For  most  problems,  we  may  obtain  a  better  physical  interpretation  of  the  results  in  the 
Fourier  transform  domain  (|)  than  in  the  function  domain  (x).  But  to  avoid  time  consuming 
computations  which  may  involve  transforms  and  inverse  transforms,  the  real  computation  should 
be  carried  out  in  the  function  domain  (x)  only.  For  the  reproducing  conditions,  Eq.(6.3.8)  is  the 
link  between  the  function  domain  (x)  and  the  Fourier  transform  domain  (^.  The  question  of  how 
to  design  a  window  function  in  the  Fourier  transform  domain  and  perform  the  computation  m 

the  function  domain  (x)  is  the  topic  of  the  next  section. 


6.4  Design  of  a  Reproducing  Kernel  Window 

Modified  Window  Function  in  the  Fourier  Transform  Domain 

A  given  window  function  may  not  satisfy  the  reproducing  conditions  and  we  have  tew 
clues  to  find  a  cure  in  the  function  domain  (x).  However,  if  we  turn  our  attention  to  th^Founer 
transform  domain  (|),  we  can  realize  that  we  can  define  a  modified  window  function,  <p^{^)  or 
0[a^),  as  a  linear  combination  of  the  original  window  function,  ^a{^)  ^(a^),  and  its  first  and 

second  derivatives,  etc.  That  is; 


II 

k=0 

or 

n 

0{a^)  =  aQ^{a^)  +  oCia^'{a^)+---+a„a''0^'^\a^)  =  ^^a,^a 

k=0 


(6.4.1a) 


(6.4.1b) 


Now  the  process  to  satisfy  the  reproducing  conditions  for  the  modified  window  function  is 
reduced  to  determining  the  proper  coefficients  <4,  i.e. 
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’  0(0) 

0'(O) 

= 

_J'"’(0). 

0(0)  fl0'(O) 


fao' 

'f 

0 

0 

_0_ 

Therefore 


'aol  r  0(0)  a0'(O)  ■••  fl>^"^(0)  ^ 

«!  0'(O)  a0"(O)  a>^"^^H0)  0 

;  =  .  ;  0 

«J  a0(""^HO)  -  LO 


(6.4.2) 


(6.4.3) 


Modified  Window  in  the  Function  Domain 

The  modified  window  function  can  be  derived  in  the  function  domain  without  applying  the  inverse 
Fourier  transform.  Since  the  reproducing  conditions  apply  at  <^=0.  an  equivalent  expression  can  be 
derived  in  the  function  domain  (x)  by  applying  the  moment  equation  defined  in  Eq.(6.3.8),  such 

that 


<  n 

0(0)  =  ^afca"0^*HO)  a^r’^  j(x  -  yf  0«(x  -  y)dy 


=  J  ^a„r\x-yf<Pa{^-y)dy  =  j(l>a{x-y)dy 


-oclk^O 


where  0(*)  is  the  modified  window  function  defined  in  the  function  domain  (x) 

n 

0a(^  -  >’)  =  X  ~ 

^=0 

In  general,  this  modified  window  function  can  be  expanded  in  the  form  of 

^^{x;x-y)  =  Po{a,x)(l>^{x-y)  +  Pi{a,x){x-y)(p^{x-y)+- 


(6.4.4) 


(6.4.5) 


+P„{a,.x){x  -  y)"  (pai^  -  >') 


6S 


n 

=  {a,x)ix  -  yf  (pa{x  -  y) 


where  the  coefficients are  associated  with  the  coefficients  «<,.  in  the  Founer  transform 
domain  (^). 


List  of  Modified  Window  Functions 

The  proper  choice  of  a  window  function,  which  satisfies  the  reproducing  condiuons  up  to 
the  order  of  polynomial  we  desire,  will  make  the  correction  function  equal  to  one.  A  list  ot 
modified  window  functions,  labeled  as  g(at).  with  the  capability  to  reproduce  the  U.x.x]  basis 

without  a  correction  function,  i.e.  Ca(r,A:-y)=l,  follow: 


Gaussian 

Box  function: 

Hat  function: 
Quadratic  spline: 
Cubic  spline: 


x-y 

a 


»-15| 
4 


a 


<t>a{x-y) 


(paix-y) 


1  1\  a 

'l3  5(x-y 


8  2V  a 

27  30rx-y 
17  17 1  a 


<Pa{^-y) 


<l>a{x-y) 


<t>a{x-y) 


(6.4.7a) 

(6.4.7b) 

(6.4.7c) 


(6.4.7d) 


(6.4.7e) 


The  list  of  modified  functions  for  the  [l,x,  x^-x^,xd]  basis,  <pa{x),is  also  given: 


Gaussian 


Box  function: 


15 


8  2 


a 


64 


8  V  fl 


a 


(Paix-y) 


225  _  525rx-yr 


4  1 


<Pa{^-y) 


(6.4.8a) 


(6.4.8b) 
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Hat  function: 
Quadratic  spline; 
Cubic  spline: 


10500 1 

f  .r  - 

“  11970| 

Cr  -  vA 

4' 

683  ' 

[  a  ) 

1  683  ' 

1  a  J 

144785  66675  f  ^  ~  >'  X  0,(,r-y) 

65216  8152  I  a  J  4076  V  a  )  \  ‘^ 

170^  _  42^r£::v V  iT^f ^ 
80347  80347  V  a  J  80347  V  a  J  J 


(6.4.8c 1 

(6.4.8d) 

(6.4.8e) 


Based  on  the  Fourier  analysis,  an  example  comparing  a  RKPM  window  function  with  its 
corresponding  wavelet  function  and  an  orthogonal  scaling  function  with  its  orthogonal  wavelet 
function  obtained  from  section  3  are  shown  in  Figure  6.1.  A  cubic  spline  function  is  used  in  this 

example. 
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In  the  frequency  domain,  the  orthogonal  scaling  function  f  (?)  serves  as  a  better  low-pass 
filter  than  the  modified  window  functions.  drC?')  and  However,  the  orthogonal  scaling 

function  contains  more  sidelobes  in  the  frequency  domain.  On  the  other  hand,  the  orthogonal 
wavelet  function  vr‘(?)  is  shown  as  a  band-pass  filter  in  the  high  frequency  teg.on  w.th  a 
sidelobe  With  the  same  frequency  range  as  the  orthogonal  wavelet  function,  the  modified  wavelet 
functions.  Wi)  and  also  serve  as  band-pass  filters  but  are  shifted  back  from  the  high 

frequency  range  without  any  sidelobe. 


6.5  Boundary  Effects  from  Finite  Domain 

When  applying  the  convolution  formula  in  a  finite  domain,  a  boundary  correction  technique 
is  -needed.  Since  the  convolution  formula  uses  the  values  of  the  neighboring  points  under  the 
support  of  window  function  to  reproduce  the  function  at  a  given  point,  neglecting  the  values  of  the 
function  outside  the  finite  domain  will  cause  error.  This  end-effect  can  be  fixed  by  enforcing  the 
reproducing  conditions  in  the  finite  domain.  The  correction  functions  change  shape  when  the 
window  function  approaches  the  boundary,  compensating  for  the  boundary  in  the  reproducing 
procedure.  Combining  the  correction  function  and  the  window  function  as  the  kernel  function. 
RKM  is  given  as  a  form  of  a  fmite  domain  convolution 

u^{x)  =  j  u(y)Qix-,x,y)(t>a{x  -  y)dy  =  |  M(y)0a(A:;x  -  y)dy  (6-5.1) 


where  0^(x;x-y)  =  is  the  kernel  function. 

6.6  Consequence  of  Discretization 

Fourier  Analysis 

The  reproducing  kernel  particle  method  (RKPM)  is  the  discrete  form  of  RKM.  The 
reproducing  conditions  and  the  criterion  to  derive  the  correction  function  are  different  from  the 
forms  in  the  continuous  system.  As  in  the  continuous  case,  we  will  first  investigate  these 
conditions  in  the  Fourier  transform  domain.  We  start  our  investigation  with  uniform  sampling.  For 
uniform  sampling,  choose 

Axj  =  Ax  and  Xj  =  nAx 


;€  Z  =  [  - -,-1,0,1, •••] 


(6.6.1) 


The  discrete  convolution  and  its  equivalent  expression  is  given  as 


4-00 


f  2mrLx\ 

,V  J 

[  Ax  J 

meZ 


0,(O)  +  (AP£r)o 


where  APET  stands  for  the  amplitude  and  phase  error  terms  and  it  is  in  the  form  of 


(6.6.2) 


/  2m7tx\ 

(A/>£r)„=X^»(^y  “  Z^  =  [....-2,-U.2,-.]  =  -;u/  (6.6,3) 

m€Z“ 

Note  that,  the  leading  term  of  Eq.  (6,6.2)  is  identical  to  that  of  the  continuous  form.  The  {APET)„ 
is  the  additional  error  introduced  by  the  discretization  procedure. 

The  reproducing  condition  for  linear  function  gives 

^  Xj  ■  <Pa{x  -  ~  (6.6.5) 

j  =  -oo 


where 


{APET)i  =  X 

meZ- 


/  2m7Cx\ 

\  A*  J 

Ax  ) 


(6.6.4) 


Similarly,  the  discrete  convolution  of  a  quadratic  function  (x^)  is  given  as 


^  x]  ■  0a(x  -  =  x:^0a(O)  ~ 

7  =  -oo 


+x\APET)q -ix{APET\+{APET)^ 


(6.6.6) 


where 


[APET).  =  Y, 

nieZ' 


Im/tx  \ 
A,x  ) 


By  induction,  the  reproducing  formula  for  an  nth  order  polynomial  yields 


y  x;  .  =  y  (-0* 

^'=0 

JZZ — OO  A.— W 


(6.6.7 ) 


(6.6.8) 


where 


(apet), = y 

meZ~ 


2m;D:  ^ 

"’zr  J 


(6.6.9) 


Comparing  the  continuous  and  the  discrete  forms,  we  find  that  the  differences  are  in  the  APETs. 
From  the  discrete  Fourier  analysis,  we  know  that  the  APETs  are  the  outcome  of  the  system 
discretization.  For  the  general  case,  the  APETs  decrease  as  the  dilation  parameter  increases,  but  the 
APETs  can  not  be  eliminated  from  the  reproducing  process.  Therefore,  we  recommend  a  quickly 
decaying  window  function  to  reduce  the  APETs  and  to  hopefully  reduce  the  error. 


7.  INTERPOLANT  ESTIMATION  AND  CONVERGENCE  OF 
REPRODUCING  KERNEL  PARTICLE  GALERKIN  METHODS 


7.1  Error  Estimation  of  Reproducing  Kernel  Methods 


The  Taylor’s  expansion  of  u{y)  with  respect  to  a  given  position  x  in  the  RKM  formula  gives 


\u{y)<f>a{x-y)dy 

=  u{x)^\{x  -  y)dy  -  m'(x)|(x  -  y)^a{x  -  y)dy+-  ■  ■ 
(x-y)"0a(x-y)tiy 

+  r~(x  -  yr'  0,{x  -  y)dy  +  h.o.t. 

(n  +  1)!  J-oo 


(7.1.1) 
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where  h.o.r  denoies  the  higher  order  terms  of  the  Taylor's  expansion.  The  equivalent  expression 
of  Eq.  (7.1.1)  by  the  definition  of  window  moments  is  given  as 


{n  +  l)\ 


(7.1.2) 


From  the  Fourier  analysis  of  RKM  and  the  reproducing  requirement  of  RKM,  the  moments  ot  the 
window  function  can  also  be  defined  as 


and 


mk{x)=  f  {x-yf(pa{x-y)dy  =  i''i>aH0) 

J  — oo 

mjt(.x)  =  f  (.X  -  yf  -  y)dy  =  k  =  l,  -  •,« 

J  — ©o 


(7.1.3) 


(7.1.4) 


The  equivalent  expression  of  Eq.(7. 1.2)  is  in  the  form  of 
U^‘‘{x)  =  u{x)  h.O.t. 

,  («  +  l)' 


(7.1.5) 


Therefore,  the  error  of  RKM  can  be  defined  as 
error{x)  =  u^‘‘  (x)  -  u{x) 

(n  +  1)! 

(n  +  l)!  t  J 

=  (0)1  +  h.  o.  t.  (7- 1  -6) 

(n  +  1)!  I 

Correspondingly,  the  interpolation  estimate  can  be  given  by  the  following  theorem. 


Theorem  7.1.1  Assume  u{x).  <p(.x)  €  C""'(n)n  where  Q  is  the  domain  of  micrcst. 

then 


u-u 


(7.1.7) 


where  a  is  the  dilation  parameter,  is  the  window  function,  C(a)  is  the  tth  denvative  ot  the 
associated  cotrecuon  function,  C,  U  a  constant  which  is  independent  of  the  dilation  parameter.  For 
RKM,  ||c(a)||^.  is  independent  of  a;  whereas  for  RKPM,  l|C(a)|^.  approaches  a  constant  when 

the  dUation  parameter  a  is  large  enough.  Note  that  r  is  defined  as  f  w  £(a),? ,  in  which  £(n)  ts 
a  normalization  factor  and  (l>'dQ  =  1 . 

Jq 

In  particular  for  /c  =  1  and  0,  we  have 


(7.1.8) 


and 


(7.1.9) 


Proof :  From  Eq.(4.2.17-19),  the  derivative  error  obtained  from  the  Taylor  expansion  can  be 
rewritten  as' 

Df‘u^‘{x)-D>‘uM  =  -  '^ru{x)yx-yfD^{it>a(x.x-y)C(x-,x-y)]dy  (7.1.10) 

a=n+l 

Taking  the  h'^  norm  on  both  sides  and  bounding  each  term  with  the  maximum  value  yields 


(7.1.11) 


For  instance,  can  be  bounded  by  l«l^^,«.*i  from  definition.  The  details  for  this  bounding 

procedure  can  be  found  in  Liu  et  al.  (1995e). 


7.2  Convergence  of  Reproducing  Kernel  Galerkin  Method 


Applying  the  above  interpolation  estimate  Eq.  (7.1.7),  we  examine  the  convergence  property  of 
the  reproducing  kernel  Galerkin  method.  To  clarify  the  concepts,  we  study  a  simple  example  of 
membrane  on  a  elastic  foundation .  The  strong  form  of  the  problem  is 
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-V^u  +  u  =  f  -ten 


(7.2.1) 


The  weak  or  variational  form  of  Eq.  (7.2.1)  is:  find  «  6  5 .  such  that  for  all  v€\) 

(7.2.2) 

a(v,M)  =  (v,/) 

where  u  is  the  trial  function,  v  is  the  test  function,/  is  the  force  term,  5  and  1)  is  the  collection 
of  trial  functions  and  test  functions,  respectively.  The  assembly  operators  a(v,«)  and  (v,/)  are 

defined  as  f ( Vv •  V«  +  vm)JQ  and  \vfdQ,  respectively.  Similarly,  the  reproducing  kernel 

Jq 

Galerkin  approximation  can  be  suted  as:  find  such  that  for  all  v  sV 

■a(v“,l;‘')  =  (w“,/) 


V-  are  Galerkin  approximaliona  which  belong  to  5“  and  here  5“  and  are  subsets  of 
5  and  V.  The  superscript  a  denotes  the  dilation  parameter  which  dictates  the  convergence  of  the 
reproducing  kernel  Galerkin  approximation. 

Remark:  There  are  distinct  differences  between  and  Recall  from  Eq.  (5.1.9), 

^  Mjuj ,  which  is  the  reconstruction  of  function  m  .  It  can  be  viewed  as  the  interpolation 

of  sampling  values  of  uixj)  denoted  by  uj:  while  is  the  reproducing  kernel  Galerkin 

approximauon.  which  can  be  defined  by  u‘‘='ZN,dj  and  dy  are  the  nodal  coefficients  derived 

from  Eq.(7.2.3). 

Theorem  7.2.1  Assume  u,  ^  n  the  convergence  rate  of  -  u\  is  governed 

by 


(7.2.4) 


where  C3  is  a  constant  which  is  independent  of  dilation  parameter  a. 


Proof.  Ul  e  =  denote  the  ettor  in  the  reproducing  kernel  particle  Galerkin  approximal.on. 

since  CZ  k).  we  may  take  v  =  v*’  in  the  variational  form  Eq.  (7.2.2)  which  gives 

(7  5) 

a(v"aO  =  (v",/) 

Subtracting  this  equation  from  Eq.  (7.2.3)  and  using  the  bUinearity  of  yield  the  well-known 

orthogonality  of  the  error, 

a(/,e)  =  0  ' 

To  acquire  the  desired  convergence  equation,  we  may  expand  aCv'"  +ey  +  e)  as  follows 


a(v^  +ey+e)  =  a(v^v^)  +  2a(v°,e)+a(e^ 

25  >0 


(7.2.7) 


>0 


so 


By  virtue  of  Eq.(7.2.6)  and  the  positive  definiteness  propeny  of  a(-,  ),  Eq.  (7.2.7)  implies 

a(e,«)<a(v‘* -i-e,v^ -l-e)  (7.2.8) 

Now  let  according  to  the  definition  of  the  variation  space  V.  Thus 

^  g  _  M  =  -  M ,  which  immediately  results  in 


a(e,e)<a(u^^  -u) 

2 

By  definition,  cl(e,^)  =  ||w^ anda(w’-w  m  — jjw-a  jj^j,then 


(7.2.9) 


||2 


(7.2.10) 


Thus  we  complete  the  proof  of  theorem  7.2.1, 

ln“-n!^.SC,n"l^|^.|uL..,lC(q4, 


(7.2.11) 


by  combing  Eq.(7.2a0)  with  theorem  7.1.1. 


Remarks:  ,  , 

1 .  The  energy  norm  of  the  error  is  now  bounded  and  the  error  estimation  and  convergence  c  e 

evaluated  based  oh  this  theorem.  j  r  u 

2.  The  convergence  study  for  problems  involving  essential  boundary  conditions  needs  urt  er 

investigation. 

7.3  Relationship  between  Interpolant  Estimation  and  Wavelet  Solution 

From  Eq.  (7. 1 .5),  the  projected  solutions  of  RKM  in  the  a-  scale  and  the  2a-scale  are  given  as 


(7.3.1) 

(7.3.2) 


U)  =  u(;t)  + 

The  RKM  wavelet  solution  in  2fl-scale  is  in  the  form  of 


Neglect  the  h.o.t’s,  the  wavelet  solution  is  proportional  to  the  error  term  given  in  Eq.  (7.1.6) 

(x)  =  (2"^^  - l)error(x)  (7.3.4) 

Eq  (7  3  4)  reveals  a  possibUity  that  the  wavelet  solution  can  be  an  mdex  for 

tile  computed  solution.  We  apply  this  error  estimate  m  adapuvity  which  is  described  next. 


8  MULTIRESOLUTION  FOR 

REPRODUCING  KERNEL  METHODS  AND  h-ADAPTIVITY 


8.1  Multiple  Scale  RKM 

The  RKM  is  a  modified  convolution  formulation  through  a  modified  window  function,  such  as 
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where 


and  n  is  the  order  of  polynomial  which  the  RKM  can  reproduce  exactly.  The  Fourier 
transformation  of  the  RKM  formula  yields 


{^)  =  u{^)i>a{i) 

=  +  (8.1.3) 

where  the  derivatives  of  window  function,  can  be  treated  as  different  filters  which 

represent  different  scales  of  resolution  in  the  RKM  formulation.  If  we  consider  the  resolution  limit 
of  every  individusl  filter,  an  equivalent  expression  of  the  filtered  response  can  be  rewritten  as 

where  Uk{^)  is  the  scale  of  the  response  which  matches  the  resolution  limit  of  the  given  filter 
Applying  the  relation  of  Eq.  (8.1.4)  to  the  RKM  formulation,  the  given  response  can  be 
separated  as 

The  inverse  Fourier  transformation  of  Eq.  (8,1.5)  is  given  in  the  form  of 


M^'-(x)=^o(-^)J  uo{y)0a{x  -  y)^y + 1  >^i{y){^~y)<t>a{^~y)^y 


r 

M, 

w  — oo 


(8.1.5) 
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«,(>■)(  jc-y)"  <>.(.':-.»)*■ 


where  u,(.t)  is  the  inverse  Fourier  transformalion  of  and  represents  the  (:-scale  ol  the  gtven 

response  u{x). 

From  the  Fourier  transform  analysis,  we  note  that  this  kernel  function  can  be  treated  as  a 
low-pass  filter  in  the  reconstruction  procedure.  The  multiple  scale  RKPM  is  defined  by  a  tamily  ot 
kernel  functions.  The  wavelets  corresponding  to  these  kernel  functions  are  defined  by. 

Therefore,  a  multiple  scale  decomposition  of  a  given  response  is  given  as. 


u\x)  =  Uq{x) 

=  VV',(x)-t-H'2(x)-t-M2(^) 

f 

I 

m 

=  ^VV,U)  +  M„(x) 

(=1 

where 


finest  scale 

two-level  decomposition 
three-level  decomposition 

m-level  decomposition 


u„ix)  =  l^C{2'"ao,x,y)<t>Jx-yMy)dy 

wjx)  =  ly/jx-y)uiy)dy 

with  Vjx  -y)  =  C(2'"-‘ao,x,y)<^„.,(x -  y)  -  C{2'"a,,x,y)<t>Jx  -  y) 


(8.1.8) 


(8.1.9) 

(8.1.10) 

(8.1.11) 


The  multiple  scale  RKPM  starts  from  the  finest  scale.  The  subsequent  levels  of  the  response  are 
obtained  by  the  decomposition  algorithm  in  Eq.(8.1.8).  Because  of  the  multiple  scale 
characteristics  of  the  time-frequency  and  space-wave  number  basis  for  representing  functions,  the 
physical  interpretation  of  the  computed  results  is  apparent. 
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8.2  Multiresolution  Analysis  and  Edge  Detection 

Edge  Location 
(wavelet  solution) 


Figure  8.1  Multiresolution  Analysis  of  Shock  Front  Structure 

Digital  edge  detection  is  an  important  technique  in  data  compression  and  image  special 
effects.  For  most  images,  edges  define  the  distinctive  features  of  the  picture.  With  the  proper 
technique,  an  image  can  be  recovered  from  its  edges  without  losing  significant  qualities.  Therefore, 
in  data  compression,  the  edges  of  images  are  used  to  optimize  the  storage  space  of  the  system  and 
to  reduce  the  translation  cost  for  network  data  exchange.  Edges  can  also  create  special  effects  in 
images.  A  blurred  image  can  be  sharpened  by  adding  distinct  edges.  On  the  other  hand,  a  picture 
can  be  smeared  by  deleting  corresponding  edges.  While  the  applications  of  edge  detection  keep 
growing,  the  mathematical  definition  of  edges  remains  rather  primitive,  with  most  researchers  still 
favoring  ad  hoc  theories  In  recent  years,  wavelet  theory  has  become  the  theory  of  choice  in  edge 

detection,  and  research  in  this  area  has  been  quite  extensive. 

In  this  section,  we  combine  multiresolution  analysis  from  RKPM  wavelet  theory  with  the 
zero-crossing  technique  to  define  the  edges  in  images.  The  procedure  is  illustrated  by  a  simple  1-D 
edge  structure  depicted  in  Fig.8.1.  Based  on  the  multiresolution  analysis,  the  edge  structure  is 
decomposed  into  two  different  scales:  the  low  scale  (the  scaling  function  solution)  and  the  high 
scale  (the  wavelet  solution).  The  low  scale  (scaUng  function)  solution  contains  the  smooth  part  of 
the  structure,  while  the  high  scale  (wavelet)  solution  provides  the  information  for  edge  detection. 
In  our  research,  the  edges  of  the  image  are  defined  as  the  zero-crossing  locations  in  the  wavelet 
solution.  In  this  1-D  case,  the  zero-crossing  position  is  located  between  the  crest  and  trough  of  the 
wavelet  like  structure  in  the  high  scale  (wavelet)  solution.  The  multiresolution  analysis  of  a  2-D 
circle  image  is  shown  in  Fig.8.2.  The  low  scale  (scaling  function)  solution  provides  an  out-of- 


Edge 

(shock  front) 


Smoothed  plot  + 

(scaling  function  solution) 


focus  or  blurred  image  of  the  original  while  the  high  scale  wavelet  solution  depicts  the  outline, 
.(edse)  of  the  circle  (image). 


Original 


Smoothed  Plot  + 


Edge  Location 


Figure  8.2  Multiresolution  Analysis  of  RKPM  in  Edge  Detection 

A  two  dimensional  edge  detection  example  is  shown  in  Fig.8.3.  In  the  first  decomposition, 
the  image  of  a  house  is  separated  into  low  and  high  (wavelet)  scales.  Applying  the  information  of 
the  wavelet  scale  and  giving  0%  and  10%  thresholds,  different  outcomes  are  observed. 


.S' 


without  threshold 


with  10%  threshold 


Figure  8.3  Edge  detection  in  image  processing 
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9.  APPLICATIONS  OF  RKPM 

9.1  Large  Deformation  Problems 

Ring  contact 

A  nonlinear  elastic  ring  contact  problem  was  studied  to  iUustrate  the  capability  ot  RKPM 
for  large  deformation  analysis.  The  material  of  the  rubber  ring  is  Mooney-Rivlin  rubber  with 
density  p  =  1.4089E-4  (sluglin^).  Other  material  properties  are  Q  =  2000. 59p5z . 
C2  =  200.361 psi  and  A  =  4E6  psi.  The  rubber  ring  impacts  the  wall  with  initial  velocity  of  1000 
(in/sec).  The  problem  statement  and  the  numerical  results  obtained  at  selected  time  instants  are 
shown  in  Fig.  9.1. 

Bar  Impact 

The  impact  of  a  two-dimensional  nonlinear  elastic  bar  was  analyzed  by  RKPM.  The 
material  for  the  bar  is  Mooney-Rivlin  rubber  with  density  p=1.4089E-4  (slug).  Other  matenal 
properties  are  Q  =  180.59p^L  =  14.67pri  and  A  =  1.47E5  psi.  A  particle  distribution  of  33 

by  9  and  timestep  of  5E-7  sec  are  used  for  computation.  Figure  9.2  shows  the  numencal  results  of 
selected  time  instants  when  the  bar  moves  toward  the  wall  with  initial  velocity  of  3000  in/sec. 

9.2  Computational  Fluid  Dynamics  (CFD) 

Under  Water  Bubble  Explosion 

Using  only  a  set  of  particles  and  Lagrangian  RKPM,  an  underwater  explosion  problem  was 
solved.  Numerical  results  describing  the  evolution  of  the  expanding  bubble  are  shown  in  Figure 
9.3.  In  the  computational  implementation,  the  cutoff  pressure  for  water  was  set  to  zero  and  the 
initial  internal  pressure  8.72E9  (Pa)  is  used.  The  evolution  of  the  expanding  bubble  is  also 
indicated  for  selected  times.  For  larger  initial  internal  pressure  case,  the  expansion  of  the  bubble 
evolved  more  rapidly  so  that  smaller  timestep  should  be  used  to  capture  the  whole  process. 

Thin  Biconvex  Airfoil  Approximation 

Two-dimensional,  inviscid,  compressible  flows  over  a  thin  biconvex  airfoil  are  examined. 
Three  different  Mach  numbers  0.5,  0.84  and  1.4  are  used  to  simulate  subsonic,  transonic  and 
supersonic  flows,  respectively.  Due  to  the  symmetry,  only  the  upper  half  of  the  airfoil  is 
considered  for  computation.  Uniform  free  stream  is  imposed  as  inflow  boundary  conditions  and  no 
values  are  prescribed  for  the  outflow  boundary  conditions.  For  a  certain  region  on  the  surface  of 
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the  airfoil.  \x\<0.5,v  = -Abxujs  specified  considering  perturbation,  where  b  is  the  relauvc 
thickness.  Press  contours  for  each  Mach  number  are  shown  in  Figure  9.4.  Aliasing  control  was 
conducted  for  comparison  and  it  can  be  seen  that  solutions  are  improved  greatly  for  transonic  a 
supersonic  flows.  Also,  2*2  multiple  scale  decomposition  for  subsonic  and  supersonic  flows  and 
3*3  multiple  scale  decomposition  for  transonic  flow  was  performed  to  separate  the  wavelet 
solution,  which  can  locate  the  high  gradient  region  more  clearly. 

Multiresolution  Analysis  and  Adaptivity 

An  appUcation  of  the  multiresolution  analysis  for  edge  detection  and  adaptivity  retmement 
in  computational  fluid  dynamics  is  iliustrated  in  Figure  9.5.  In  this  multiresolution  analysis,  the 
numerical  solution  is  examined  through  error  estimation  and  convergence  studies.  In  the  multilevel 
refinements,  only  a  set  of  nodes  are  required  for  adaptivity  with  no  apprehension  ot  compatibility 
problems.  In  addition,  the  physical  interpretation  of  the  numerical  solutions  can  be  synthesized  by 

multiple  scale  representation. 

9.3  Structural  Acoustics 
Plane  Wave  Scattering  Problems 

An  incident  wave  <^9  travels  from  left  to  right  along  the  x-axis  (6=0)  with  a  rigid  circle 
cylinder  appearing  on  its  path.  The  geometry  and  problem  setup  are  shown  in  Fig.9.6a. 

The  solution  domain  is  the  shaded  area  between  the  inner  boundary,  representing  the  cylinder,  and 
the  outer  artificial  boundary,  simulating  the  infinite  domain.  The  first  order  DtN  boundary 

condition. 


is  used  in  the  computation.  And  the  rigid  surface  reflection  has  the  form. 


<P.n  +  =0  on  Fh 


The  series  solution  is  given  as: 
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t,  =  -2'£ 


■^■j'?LffW(/cr)cos(iie) 


Figure  9.6a  Rigid  infinite  circular  cylinder 

The  numerical  results  are  shown  in  Fig.9.6b.  Note  that  the  edge  detection  technique  is  also 
used  to  conduct  hp-adaptivity  in  this  example. 

Multiple  Scale  Analysis 


The  flowchart  of  multiple  scale  analysis  for  structural  acoustics  is  outlined  in  Fig.9.7.  By 
introducing  an  adequate  filter,  original  solutions  can  be  separated  into  low  scale  and  high  scale 
solutions  (also  known  as  wavelet  solutions),  which  will  provide  better  physical  interpretation  of 
the  results.  Applying  the  zero-crossing  technique  described  in  Section  8.2  on  the  wavelet  solution, 
the  edges  of  the  response  can  be  determined.  Since  the  edges  represent  the  high  gradient  region  of 
the  solution,  they  automatically  provide  the  locations  for  local  refinement.  To  enhance  the 
resolution  in  high  gradient  region,  extra  particles  are  inserted  on  the  edges.  The  corresponding 
window  functions  are  adjusted  to  match  the  new  distribution  of  the  particles.  This  procedure  is 
similar  to  the  classical  hp-adaptive  refinement  and  the  behavior  of  the  solutions  is  also  similar.  This 
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kind  of  refinement  can  continue  widtout  encounmring  any  difficulties  until  the  resolution  lint,.  ,s 
reached,  i.e.  desired  solution  is  obtained.  A  numerical  example  was  studied  using  this  procedure 
and  the  results  are  given  in  Fig.9.7.  For  die  starting  level,  684  nodal  points  were  used.  Four  eve  .s 
of  refinements  were  performed  to  acquire  the  satisfactory  results  using  2407  nodal  points. 

Green's  Function  for  a  Rectangular  Domain 

Consider  the  problem  of  fmding  the  Green’s  function  with  homogeneous  Dmchlet 
boundary  condition  within  a  rectangular  domain  (see  figure  9.8a). 


V2«+>c^«  =  / 

where  ir^  =  (m/cf  and  /  =  5(x  -  x„ -  yo)  ( =  .Vo  *  0.5  )wi,h  boundary  conditions 
w(0,  y)  =  u{L,  y)  =  u(x,0)  =  u{x,  L)  =  0 


The  series  solution  to  this  problem  is  given  as; 


«  «  sm 

m- 


mJtXr 


sin 


mitx 

~~r 


Note  that  there  exist  natural  frequencies  whenever 


2 
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A  uniform  25x25  particles  with  mesh  size  Ax=Ay=l/24  is  used  for  RKPM  solution.  It  is 
equivalent  a  24x24  linear  FEM  mesh,  a  12x12  quadratic  FEM  mesh,  or  a  8x8  cubic  FEM  mesh. 
The  cubic  spline  function  is  used  as  the  window  function  in  the  RKPM  with  dE/E=0.1%  as  the 
aliasing  control  parameter  [Liu  and  Chen  (1995)].  The  resolution  limit  is  JC(wave  number)=75.40 
or  ;r/  (1/  24).  The  La-norm  of  the  series  and  numerical  solutions  with  the  wave  numbers  swept 
from  0  to  80  are  shown  in  Fig.  9.8b.  When  the  wave  number  is  approached  to  the  natural 
frequency  of  the  system,  we  can  observe  a  peak  in  the  La-norm  plot.  If  the  discrete  system  is  a 
good  approximation  of  the  continuous  system,  the  peaks  from  two  systems  will  match  point  to 
point  in  the  frequency  (wave  number)  domain.  In  the  low  wave  number  region,  those  peaks  match 
quiet  weU  for  the  most  of  the  cases,  but  begin  to  fail  to  match  with  the  continuous  system  in  higher 
wave  number  region.  The  position  where  the  peaks  begin  to  separate  will  determine  the  resolvable 
scale  of  the  given  system  (the  combination  effect  of  discretization  and  interpolation).  As  depicted  in 
the  plots,  the  linear  finite  element  fails  quite  early,  and  then  follows  the  quadratic  FE.  The  cubic  FE 
does  not  survive  long  enough  either.  The  RKPM  can  extend  its  resolvable  wave  number  up  to  60, 
almost  80%  (60/75.4)  of  the  resolution  limit  In  other  words,  we  can  use  only  3.5  points  to  catch  a 
cos-sin  wave  exactly  by  RKPM  (the  resolution  limit  is  3.14  points!). 


0.0044  sec 


iiiiiiiiiiiiiiiiiiiiiiiiiii  5.oe.5.« 


l.OE-4  sec 


\.5E-Asec 


2.0E-4  sec 


2.5E-4  sec 


Figure  9.2  Bar  impact  for  large  deformation  problem 
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Figure  9.3  Under  water  bubble  explosion 


Mach,  number  =  0.5  (Subsonic) 


high  wave  number  band  (wavelet)  solutions 

RKPM  Solutions  Without  Aliasing  Control  RKPM  Solutions  With  Aliasing  Control  depicting  the  shock  (or  high  gradient)  locations 


scalingpnction 


Figure  9.5  Mullirc.soliilion  analy.si.s  and  adaplivily  l\)i  CI  D 


Figure  9.6b  Plane  wave  scattering  problems 
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Figure  9.7  Multiple  scale  analysis  for  structural  acoustics 
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CONCLUSIONS 

The  reproducing  kernel  particle  method  (RKPM)  was  reformulated  in  the  framework  of 
wavelet  theories  which  include  such  concepts  as  reproducing  conditions,  discrete  convolutions, 
and  multiresolution  analysis.  It  is  shown  that  for  a  sufficiently  smooth  function,  the  interpolant 
expansion  in  terms  of  sampled  values  will  converge  to  the  original  function  in  Sobolev  norm.  With 
this  construction  of  interpolation  error,  the  convergence  rate  can  be  measured  in  terms  a  new 
control  variable,  the  dilation  parameter  of  the  window  function. 

By  studying  the  convolution  of  spline  family  with  polynomials,  we  are  able  to  construct  a 
new  family  of  compactly  supported  scaling  functions  which  can  reconstruct  polynomials  of 
arbitrary  order  exactly.  It  was  first  shown  that  the  interpolant  error  estimate  is  proportional  to  the 
wavelet  solution.  We  have  also  implemented  the  first  application  of  wavelet  solution  and  edge 
detection  in  designing  hp  adaptive  algorithm. 

RKPM  has  been  applied  to  such  areas  as  structural  acoustics,  elastic -plastic  deformation, 
computational  fluid  dynamics.  RKPM  is  shown  to  provide  an  accurate  mesh  free  algorithm  and 
possess  superior  convergence.  RKPM  is  envisioned  as  a  new  meshless  method  with  completely 
different  mentality,  which  will  build  a  bridge  between  the  traditional  interpolation  method  and  the 
spectral  method.  This  is  very  desirable  for  the  next  generation  of  meshless  methods. 
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appendix  a  3D  SHAPE  FUNCTION  AND  DERIVATIVE 


Al.  Shape  Function 

For  a  linear  basis  in  3D,  the  reproducing  equation  is  (c.f.  Eq.  5.3.3): 


NP 

(X)  =  •^u(xj)C,(x-.x  -Xj)i>.{x- Xj )hVj 


(Al.l) 


y=i 


.  ,  ,  t X,  X,  \th  (7  -  7  ^  (b  (b  X  <t>n  are  cubic  spline  functions, 

where  (x  -Xj)  =  <ba,  (x-Xj  )(t>a^  (y  -  yj  )<pa,  U  Zj),  (Pa,  X  (Pa^  X 

and  AV  is  the  nodal  volume.  The  correction  function  Caix-,x-Xj)  is  defmed  as: 


Ca(x;x-xj)  =  boix)  +  bi{x){x  -  xj)  +  b2{x)(y  -  yj)  +  hix)(z  -  zj) 


The  correction 


function  coefficients  Hx)  are  computed  by  Eq.  (5.3.6): 


b(x)  =  M  \a,x)PiO) 


(A1.2) 


(A1.3) 


where 


F(0)  =  [1,0,0,  of 


M{a,x)  = 


niQQQiayX)  mx^QQ{a,x)  ;^io(fl,Jf) 
mi(X)(«,Jc)  /^oo(fl,Jc)  miio(®»'^) 

;noio(a,jc)  mi  10 (<*>*)  ^20 

niQQiiayX)  mioi(a,x)  /^ii(a,x)  moo2(<*>*) 


and 


NP 


(A1.4a) 


(A  1.4b) 


(A  1.4c) 


=  2,  ix-xj)‘{y-yj)\z  -  Zjf  0,(x-xj)AVj 

)=1 

The  shape  function  is  then  given  as: 

Nj{x-,x-Xj)  =  Ca{x\x-Xj)4>„{x-Xj)AVj 

=  [bo{x)  +  biix){x-Xj)  +  b2(x){y-yj)  +  bi  {x)iz  -  Zy)] 

(A1.5) 


A2.  First  Derivative  of  Shape  Function 

The  first  derivative  of  shape  function  is  given  as: 


N,  {x:x-x.)  =  \Ca,{x;x-Xj)0aix-Xj)  +  Ca(x-,x-Xj)^a.Ax-Xj)]AVj 

J,X  J  I 


{A2A) 


or  in  matrix  form: 


Ca  xix\X  -  X  j)0a{x  -  Xj)+  Ca{x;X  -  Xj)^a^^{X  -  Xj) 

Ca.y{x;x  -  Xj)^a(^  ~  ~  ~ 

Ca  AX\X  -  X j)<!>a(x  -  Xj)  +  Ca(x;X  -Xj)<Pa  ^{X  -Xj) ^ 


(A2.2) 


where 


C^^ix;x-Xj)  =  bo_y{x)-^b,^y{x)(x-xj)  +  b2^y{x){y-yj)  +  t^^^^^^^^ 

Cajx-,x  -Xj)  =  b^  Jx)  +  foi  ,(x)(x  -Xj)  +  b2  jx){y  -  Vy )  +  b^Jx)iz  -  Zj)  +  hix) 

The  first  derivative  of  the  correction  function  coefficients  are  obtained  as  follows. 

; 

b^(x)  =  -M~\a,x)M^^ia,x)b{x)- 


C  C  J 
fu 

where 


-n  o 

j  Ci  '. 


'»Um^(a,x')  mioo,*(«>Jc)  m^i,,ia,x) 

'"too, miiQj^ia,x)  m^oi^ittyx) 
mQQi^(a,x)  mioij,(a,Jf)  ^Uj,(a,x)  /noo2^(«>-^)_ 


(A2.3a) 


(A2.3b) 


(A2.3c) 


(A2.4) 


(A2.5a) 


(A2.5b) 


